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Abstract. Several issues related to the gravitational clustering of collision- 
less dark matter in an expanding universe is discussed. The discussion is 
pedagogical but the emphasis is on semianalytic methods and open ques- 
tions — rather than on well established results. 



1. Mathematical description of gravitational clustering 

The gravitational clustering of a system of collisionless point particles in an 
expanding universe poses several challenging theoretical questions. Though 
the problem can be tackled in a 'practical' manner using high resolution 
numerical simulations, such an approach hides the physical principles which 
govern the behaviour of the system. To understand the physics, it is neces- 
sary that we attack the problem from several directions using analytic and 
semianalytic methods. These lectures will describe such attempts and will 
emphasise the semianalytic approach and outstanding issues, rather than 
more well established results. In the same spirit, I have concentrated on the 
study of dark matter and have not discussed the baryonic physics. 

The standard paradigm for the description of the observed universe pro- 
ceeds in two steps: We model the universe as made of a uniform smooth 
background with inhomogeneities like galaxies etc. superimposed on it. 
When the distribution of matter is averaged over very large scales (say, 
over 200 h~^Mpc) the universe is expected to be described, reasonably 
accurately, by the Friedmann model. The problem then reduces to under- 
standing the formation of small scale structures in this, specified Friedman 
background. If we now further assume that, at some time in the past, there 
were small deviations from homogeneity in the universe then these devia- 
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tions can grow due to gravitational instability over a period of time and, 
eventually, form galaxies, clusters etc. 

The study of structure formation therefore reduces to the study of the 
growth of inhomogeneities in an otherwise smooth universe. This — in turn 
— can be divided into two parts: As long as these inhomogeneities are small, 
their growth can be studied by the linear perturbation around a background 
Friedmann universe. Once the deviations from the smooth universe become 
large, linear theory fails and we have to use other techniques to understand 
the nonlinear evolution. [More details regarding structure formation can be 
found e.g. in Padmanabhan, 1993; 1996] 

It should be noted that this approach assumes the existence of small 
inhomogeneities at some initial time. To be considered complete, the cos- 
mological model should also produce these initial inhomogeneities by some 
viable physical mechanism. We shall not discuss such mechanisms in these 
lectures and will merely postulate their existence. There is also a tacit 
assumption that averaging the matter density and solving the Einstein's 
equations with the smooth density distribution, will lead to results com- 
parable to those obtained by averaging the exact solution obtained with 
inhomogeneities. Since the latter is not known with any degree of confi- 
dence for a realistic universe there is no straightforward way of checking 
this assumption theoretically. [It is actually possible to provide counter ex- 
amples to this conjecture in specific contexts; see Padmanabhan, 1987] If 
this assumption is wrong there could be an effective correction term to the 
source distribution on the right hand side of Einstein's equation arising 
from the averaging of the energy density distribution. It is assumed that 
no such correction exists and the universe at large can be indeed described 
by a Friedmann model. 

The above paradigm motivates us to study the growth of perturbations 
around the Friedmann model. Consider a perturbation of the metric gaf}{x) 
and the stress-tensor Tap into the form {gap+Sigap) and {Tai3+ 8X^13), where 
the set {gap-iTafj) corresponds to the smooth background universe, while 
the set (Sgaf^jSTais) denotes the perturbation. Assuming the latter to be 
'small' in some suitable manner, we can linearize Einstein's equations to 
obtain a second-order differential equation of the form 

^9ap)89al3 = 8Ta0 (1) 

where £ is a linear differential operator depending on the background space- 
time. Since this is a linear equation, it is convenient to expand the solution 
in terms of some appropriate mode functions. For the sake of simplicity, let 
us consider the spatially flat [Vl = 1) universe. The mode functions could 
then be taken as plane waves and by Fourier transforming the spatial vari- 
ables we can obtain a set of separate equations -C(k)'55'(k) = '^^(k) for each 
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mode, labeled by a wave vector k. Here is a linear second order differ- 
ential operator in time. Solving this set of ordinary differential equations, 
with given initial conditions, we can determine the evolution of each mode 
separately. [Similar procedure, of course, works for the case with 0,^1. 
In this case, the mode functions will be more complicated than the plane 
waves; but, with a suitable choice of orthonormal functions, we can obtain 
a similar set of equations] . This solves the problem of linear gravitational 
clustering completely. 

There is, however, one major conceptual difficulty in interpreting the 
results of this program. In general relativity, the form (and numerical value) 
of the metric coefficients (or the stress-tensor components Tap) can 
be changed by a relabeling of coordinates — a;"'. By such a trivial 
change we can make a small dT^p large or even generate a component 
which was originally absent. Thus the perturbations may grow at different 
rates — or even decay — when we relabel coordinates. It is necessary to 
tackle this ambiguity before we can meaningfully talk about the growth of 
inhomogeneit ies . 

There are two different approaches to handling such difficulties in gen- 
eral relativity. The first method is to resolve the problem by force: We 
may choose a particular coordinate system and compute everything in that 
coordinate system. If the coordinate system is physically well motivated, 
then the quantities computed in that system can be interpreted easily; 
for example, we will treat 5Tq to be the perturbed mass (energy) density 
even though it is coordinate dependent. The difficulty with this method is 
that one cannot fix the gauge completely by simple physical arguments; the 
residual gauge ambiguities do create some problems. 

The second approach is to construct quantities — linear combinations of 

various perturbed physical variables — which are scalars under coordinate 
transformations, [see eg. the contribution by Brandenberger to this vol- 
ume and references cited therein] Einstein's equations are then rewritten 
as equations for these gauge invariant quantities. This approach, of course, 
is manifestly gauge invariant from start to finish. However, it is more com- 
plicated than the first one; besides, the gauge invariant objects do not, in 
general, possess any simple physical interpretation. In these lectures, we 
shall be mainly concerned with the first approach. 

Since the gauge ambiguity is a purely general relativistic effect, it is 
necessary to determine when such effects are significant. The effects due to 
the curvature of space-time will be important at length scales bigger than 
(or comparable to) the Hubble radius, defined as ^^(t) = (d/a)~^. Writing 
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the Friedmann equation as 

a2 



(2) 



where ^r, Unr, and represent the density parameters for relativistic 
matter (with pR = [1/3)pr; pr oc a~^), non relativistic matter with pnr = 
0; Pnr oc a~^), cosmological constant {pv = —pv'-, pv = constant) and total 
energy density = Q.r + ^nr + ^v)-, respectively, it follows that 



dH{z) = H^^ \p.R{l + zf + J^nr(1 + zf + (1 - + zf + J^y 
This has the limiting forms 



-1/2 

(3) 



1 H^'n^K^l + ^)-3/2 » z » Zeurv; = 0) 

during radiation dominated and matter dominated epochs where 

(1 + Zeg) = ^^; {l + Zcurv)^T^ 1 (5) 

(The universe is radiation dominated for z S> Zgq and makes the transition 
to matter dominated phase at z ~ z^q. It becomes 'curvature dominated' 
sometime in the past, for z ^ Zcuw, if ^NR < 0.5. We have set ily = 
for simplicity). The physical wave length Aq, characterizing a perturbation 
of size Aq today, will evolve as A(z) = Ao(l + z)^^ t. Since dn increases 
faster with redshift, (as (1 + z)^^/^ in matter dominated phase and as 
(1 + z)~'^ in the radiation dominated phase) X{z) > dn^z) at sufficiently 
large redshifts. For a given Aq we can assign a particular redshift Zgnter at 
which A(zcntcr) = '^//(-^cntcr)- For z > Zcntcr, the proper wavelength is bigger 
than the Hubble radius and general relativistic effects are important; while 
for z < Zenter wc have \ < dn and one can ignore the effects of general 
relativity. It is conventional to say that the scale Aq "enters the Hubble 
radius" at the epoch Center- 

The exact relation between Aq and Center differs in the case of radiation 
dominated and matter dominated phases since dniz) has different scalings 
in these two cases. Using equation (4) it is easy to verify that: (i) A scale 

Aeq ^ (^j in](^/nNR) = 14Mpc(J^NR/i')-' (6) 
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enters the Hubble radius at z = z^q. (ii) Scales with A > Aeq enter the 
Hubble radius in the matter dominated epoch with 



(iii) Scales with A < Aeq enter the Hubble radius in the radiation dominated 
epoch with 



One can characterize the wavelength Aq of the perturbation more mean- 
ingfully as follows: As the universe expands, the wavelength A grows as 
X{t) = Xo[a{t)/ao] and the density of non-relativistic matter decreases as 
p(t) = po[ao / a{t)]^ . Hence the mass of nonrelativistic matter, M(Ao) con- 
tained inside a sphere of radius (A/2) remains constant at: 




This relation shows that a comoving scale Aq ~ 1 Mpc contains a typical 
galaxy mass and Aq ~ 10 Mpc contains a typical cluster mass. Prom (8), we 
see that all these — astrophysically interesting — scales enter the Hubble 
radius in radiation dominated epoch. 



This feature suggests the following strategy for studying the gravita- 
tional clustering. At 2; ^ Center (for any given Aq), the perturbations need 
to be studied using general relativistic, linear perturbation theory. For 
z <C ^^cntcD general relativistic effects are ignorablc and the problem of 
gravitational clustering can be studied using newtonian gravity in proper 
coordinates. Observations indicate that the perturbations are only of the 
order of (10"^ - IQ-^) at z — ^cntcr for all Aq. Hence the nonlinear epochs 
of gravitational clustering occur only in the regime of newtonian gravity. In 
fact the only role of general relativity in this formalism is to evolve the ini- 
tial perturbations upto z ^ -Zcnter; after which newtonian gravity can take 
over. Also note that, in the nonrelativistic regime {z ^ Zcntcr : ^ dn), 
there exists a natural choice of coordinates in which newtonian gravity is 
applicable. Hence, all the physical quantities can be unambiguously defined 
in this context. 

2. Linear growth in the general relativistic regime 

Let us start by analysing the growth of the perturbations when the proper 
wavelength of the mode is larger than the Hubble radius. Since A S> dn we 
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cannot use newtonian perturbation theory. Nevertheless, it is easy to deter- 
mine the evolution of the density perturbation by the following argument. 

Consider a spherical region of radius containing energy den- 

sity pi = Pb + Sp, embedded in a A: = Priedmann universe of density pb- 
It follows from spherical symmetry that the inner region is not affected by 
the matter outside; hence the inner region evolves as a A; 7^ Priedmann 
universe. Therefore, we can write, for the two regions: 

H^ = ^Pb, H^ + ^ = ^{pb + Sp). (10) 

The change of density from pb to pb + Sp is accommodated by adding a 
spatial curvature term (k/a^). If this condition is to be maintained at all 
times, we must have 

-^Sp = -^, (11) 

or 

^ = ^ (12) 

Pb SnGipba^y ^ ' 

If {6p/pb) is small, a{t) in the right hand side will only differ slightly from 
the expansion factor of the unperturbed universe. This allows one to deter- 
mine how {6p/pb) scales with a for A > dn- Since pb oc in the radiation 
dominated phase {t < teq) and pb oc in the matter dominated phase 
{t > teq) we get 

pj \a (for t> teq). ^ ^ 

Thus, the amplitude of the mode with A > cIh always grows; as in 
the radiation dominated phase and as a in the matter dominated phase. 
Since no microscopic processes can operate at scales bigger than dn all 
components of density (dark matter, baryons, photons), grow in the same 
manner, as (5 oc (p^a^)"^ when A > dff. 

A more formal way of obtaining this result is as follows: We first recall 
that there is an exact equation in general relativity connecting the geodesic 
acceleration g with the density and pressure: 

V • g = -inGip + 3p) (14) 

Perturbing this equation, in a medium with the equation of state p = wp, 
we get 

Vr • [Sg] = -AttG {5p + 35p) = -47rG/9fe (1 + Su-) 5 = a'^Vx • [5g] (15) 

where 5 = (5p/ p) is the density contrast. Let us produce a (5g by introducing 
a perturbation of the proper coordinate r = a(t)x to the form r -|- 1 = 
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a(i)x[l + e] such that 1 = axe. The corresponding perturbed acceleration is 
given by Sg = x[oe + 2ae]. Taking the divergence of this Sg with respect to 
X we get 

Vx • [5g] = 3 [ae + 2de] = -4TrGpba{l + 3w)5 (16) 
This perturbation also changes the proper volume by an amount 

{SV/V) = {Sl/r) = 3e (17) 

If we now consider a metric perturbation of the form gik — > gik + hik, the 
proper volume changes due to the change in by the amount 

(SV/V) = -ih/2) (18) 

where h is the trace of hik- Comparison of the expressions for (dV/V) 
suggests that, as far as the dynamics is concerned, the equation satisfied by 
3e and that satisfied by —{h/2) will be identical. Substituting e = {—h/6) 
in equation (16), we get 

h + 2(^^ h = 87rGp6(l + 3w)S (19) 

(A more formal approach — using full machinery of general relativity — 
leads to the same equation.) We next note that 6 and h can be related 
through conservation of mass. Prom the equation d{pV) = —pdV we obtain 

S=^ = -il + w)^ = -3{l + w)e (20) 

giving 

6 = -3e{l + w) = +{l + w)^ (21) 
Combining (19) and (21) we find the equation satisfied by 6 to be 

5 + 2'^5 = AT^Gpb{l + w){l + 3w)6. (22) 

This is the equation satisfied by the density contrast in a medium with 
equation of state p = wp. 

To solve this equation, we need the background solution which deter- 
mines a{t) and Pb{t). When the background matter is described by the 
equation of state p = wp, the background density evolves as pi, oc a~^^^^'^\ 
In that case, Friedmann equation (with Q = 1) leads to 
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5g<xf- 5dOc-] n = - \J^ (25) 



provided w 7^ —1. When w = —1, a{t) oc exp(//i) with a constant We 
will consider w ^ —1 case first. Substituting the solution for a{t) and pb{t) 
into (22) we get 

c, 4 ^_2 (l + 3u;) <5 

3(1 + u;) t 3 (1 + u;) ^ ^ 

This equation is homogeneous in t and hence admits power law solutions. 
Using an ansatz S oc t", and solving the quadratic equation for n, we find 
the two linearly independent solutions {dg,6ci) to be 

1 _ _2{l + 3w) 
t' ^ ~ 3 (l + w) 

In the case of li; = — 1, a(i) oc exp (fit) and the equation for 6 reduces to 

S + 2X6 = 0. (26) 

This has the solution 5g oc exp(— 2/it) oc a"^. All the above solutions can 
be expressed in a unified manner. By direct substitution it can be verified 
that 6g in all the above cases can be expressed as 

<5g oc (27) 

which is exactly the result obtained originally in (12). This allows us to 
evolve the perturbation from an initial epoch till z = -Center, after which 
newtonian theory can take over. 



3. Gravitational clustering in Newtonian theory 

Once the mode enters the Hubble radius, dark matter perturbations can be 
treated by newtonian theory of gravitational clustering. Though Sx<^l at 
z ^ ;Zcntcrj we shall develop the full formalism of newtonian gravity at one 
go rather than do the linear perturbation theory separately. 

In any region small compared to dn one can set up an unambiguous co- 
ordinate system in which the proper coordinate of a particle r{t) = a(t)x(t) 
satisfies the newtonian equation = — Vr$ where ^ is the gravitational 
potential. Expanding r and writing $ = $frw + where ^I'frw is due to 
the smooth component and (f) is due to the perturbations, we get 

ax + 2dx + ax = - Vr^FRW - Vr<?!' = - Vr^FRV^ - a~-^Vx<?!' (28) 

The first terms on both sides of the equation (ax and — Vr^FRw) should 
match since they refer to the global expansion of the background FRW 
universe. Equating them individually gives the results 

d. 1„, , 1 a n 27rG , „ x 9 

x + 2-x = -— Va;<^ ; $FRW = -TT-r = ^ P + 3p r2 29 

a a^ 2 a 6 
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where (p is generated by the perturbed, newtonian, mass density through 

Vlcf) = 47rGa^ {5p) = ATrGpbO^S. (30) 

If Xi(t) is the trajectory of the i — th particle, then equations for newtonian 
gravitational clustering can be summarized as 

±i + —±i = -\v^cP; Vl<P = A7rGa^pbS (31) 
a 

where ph is the smooth background density of matter. We stress that, in the 
non-relativistic limit, the perturbed potential (p satisfies the usual Poisson 
equation. 

Usually one is interested in the evolution of the density contrast S {t, x) 

rather than in the trajectories. Since the density contrast can be expressed 
in terms of the trajectories of the particles, it should be possible to write 
down a differential equation for (5(t,x) based on the equations for the tra- 
jectories Xj(t) derived above. It is, however, somewhat easier to write down 
an equation for (5k(t) which is the spatial fourier transform of 6{t, x). To do 
this, we begin with the fact that the density p(x, t) due to a set of point 
particles, each of mass m, is given by 

^^""^ = ^ ^ '^^^'^ " ^^^^ 

where Xi(i) is the trajectory of the ith particle. To verify the normal- 
ization, we can calculate the average of /o(x, t) over a large volume V. We 
get 

where N is the total number of particles inside the volume V and M = Nm 
is the mass contributed by them. Clearly pi, oc a~^, as it should. The density 
contrast (5(x, t) is related to /^(x, t) by 

l + <^(x,i)^^i^ = ^^5^[x-x,(i)]= /dVD[x-xr(i,q)]. (34) 

In arriving at the last equality we have taken the continuum limit by re- 
placing: (i) Xj(t) by XT(t, q) where the initial position q of a particle lables 
it; and (ii) (V/N) by d^q since both represent volume per particle. Fourier 
transforming both sides we get 



Sk{t) = / (f:ic^^''S{x,t) = / (fci exp[-ik.xr(t,q)] - {2Trf Sd{^) (35) 
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Differentiating this expression, and using the equation of motion (31) for 
the trajectories give, after straightforward algebra, the equation: 



(5k + 2-(5k = 47rGp6(5k + ^k - 
a 



with 



^k = 4:710 pb J ^^;^(5k"5k-k 



(27r)2 



k.k' 



-Bk = y" c^'^q (k.xr)^ cxp [-ik.XTit, q)] . 



(36) 



(37) 



(38) 



This equation is exact but involves xriti q) on the right hand side and hence 
cannot be considered as closed, [see, eg. Peebles, 1980; the expression for 
Ak is usually given in symmetrised form in k' and (k— k') in the literature]. 

The structure of (36) and (38) can be simplified if we use the perturbed 
gravitational potential (in Fourier space) 0k related to 6^ by 



3Hi 



and write the integrand for Ai^ in the symmetrised form as 



(39) 



k'Ok-k' 



k.k' 



k'^ 



-^k'<^k-k' 



1 /Si 



k.k' ^ k.(k-k') 



/t'2 |k-k'|2 
f^k-k' 



2 Kk'^J V Ik-k'|2 



1 / 2a 



2 \3H^ 



9k'^^k-k' 



(k-k')2k.k' + /c'2 (fc2-k.k') 
fe2(k.k' + fc'2) -2(k.k')2 



(40) 



In terms of 0k, equation (36) becomes, for a O = 1 universe. 



A ; 1 
&k + 4-0k = -7r2 
a 2a^ 



!'k'<Pk-k' 



k'.(k + k')-2 



k.k 



+ 



3^\ f<l^^(^\\^. 
2 J a \ k 



(41) 



where x = xr(t, q). We shall see later how this helps one to understand 
power transfer in gravitational clustering. 
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If the density contrasts are small and linear perturbation theory is to be 
valid, we should be able to ignore the terms and in (36). Thus the 
liner perturbation theory in newtonian limit is governed by the equation 



+ 2-jk 
a 



(42) 



From the structure of equation (36) it is clear that we will obtain the linear 
equation if ^k ^ 4:7rGpb5k and Bk ^ 47rGpft(^k- A necessary condition for 
this dk 1 but this is not a sufficient condition — a fact often ignored or 
incorrectly treated in literature. For example, if (5k ^ for certain range of 
k at t = to (but is nonzero elsewhere) then A^^ ^ AirGphdi^ and the growth 
of perturbations around k will be entirely determined by nonlinear effects. 
We will discuss this feature in detail later on. For the present, we shall 
assume that ^k and Sk are ignorable and study the resulting system. 

4. Linear perturbations in the Newtonian limit 

At z ^ Center) the perturbation can be treated as linear (p ^ 1) and new- 
tonian (A <C d//). In this case, the equations are 



Sk + 2-Sk = 4:TrGpDMSk 
a 



(43) 



d2 



+ 



k SttG 



{PR + PDM + Pv) 



(44) 



where poM, PR, and pv are defined in section 1. We will also assume that the 
dark matter is made of collisionless matter and is perturbed while the en- 
ergy densities of radiation and cosmological constant are left unperturbed. 
Changing the variable from t to a, the perturbation equation becomes 



2a^ 



PR + PDM + Pv - 



3k 



£6 



+ a 



f 3k \ 
2PR + 3pnM + 6py - 4 i^-^ 



-y- = '^PDmO 
da 



(45) 



Introducing the variable r = (a/flo) = and by writing pi = QiPc for 

the i^^ species, and k = — (87rG/3)pcCio(l — we can recast the equation 
in the form 



2r 



nvT^ + (1 - + ^DMT + 



5" 
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+ 



60yr^ + 4:{l-n)T^ + SnDMT + 20r 6' = SUdmS 

(46) 



where the prime denotes derivatives with respect to r. This equation is in 
a form convenient for numerical integration from r = Tenter = (1 + -Center )~^ 
to r = 1. 

The exact solution to (46) cannot be given in terms of elementary func- 
tions. It is, however, possible to obtain insight into the form of solution by 
considering different epochs separately. 

Let us first consider the epoch 1 ^ z ^ -Zenter when we can take = 
0,^ = 1, reducing (46) to 

2t (noMT + flu) 5" + CHIdmt + 2^r) 5' = SUdmS (47) 
Dividing thoughout by and changing the independent variable to 

:, = r(^)= ^ = ^ (48) 

we get 

2a;(l+a:)^ + (2 + 3x)^ = 35DM; x=^. (49) 

ax^ ax Oeq 

One solution to this equation can be written down by inspection: 

<5dm = 1 + ^x. (50) 

In other words (5dm ~ constant for a <^ Ocq (no growth in the radiation 
dominated phase) and 5du (x o for a ^ Oeq (growth proportional to a in 
the matter dominated phase). 

We now have to find the second solution. Given the first solution, the 
second solution A can be found by the Wronskian condition {Q'/Q) = 
-[(2 + 3x)/2x{l + x)] where Q = (5dm A' - (^dm^- Writing the second 
solution as A = f{x)Sj:,M{x) and substituting in this equation, we find 

f ^ _ 2 + 3x 

/' (5dm 2x(l + x)' ^ ' 

which can be integrated to give 

•^"~/x(l + 3x/2)2(l + x)V2- 
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The integral is straightforward and the second solution is 



A = /(5dm = + Y 



In 



(l + x)V^ + l 
(l + x)V2_i 



3(1 (53) 



Thus the general solution to the perturbation equation, for a mode which 
is inside the Hubble radius, is the linear superposition 6 = AS^m + BA 
with the asymptotic forms: 

X f\ AX ^^IHA^^ (A + Bln{A/x) (x <Cl) 

S,en{x) = A6j,m{x) + BA{x) = I ^3/2)Ax + (4/5)i?x(-3/2) (a; » 1). 

(54) 

This result shows that dark matter perturbations can grow only logarith- 
mically during the epoch Oenter < a < Oeq. During this phase the universe 
is dominated by radiation which is unperturbed. Hence the damping term 
due to expansion {2a/a)6 in equation (42) dominates over the gravitational 
potential term on the right hand side and restricts the growth of perturba- 
tions. In the matter dominated phase with a 3> a^q, the perturbations grow 
as a. This result, combined with that of section 2, shows that in the matter 
dominated phase all the modes (ie., modes which are inside or outside the 
Hubble radius) grow in proportion to the expansion factor. 

Combining the above result with that of section 2, we can determine 
the evolution of density perturbations in dark matter during all relevant 
epochs. The general solution after the mode has entered the Hubble radius 
is given by (54). The constants A and B in this solution have to be fixed by 
matching this solution to the growing solution, which was valid when the 
mode was bigger than the Hubble radius. Since the latter solution is given 
by 6{x) = x^ in the radiation dominated phase, the matching conditions 
become 

2Xenter = [A5'j,m{x) + B A' {x)] ^^^^^^^^ . 

(55) 

This determines the constants A and B in terms of Xcnter = ( Center /oieq) 
which, in turn, depends on the wavelength of the mode through Oenter- 

As an example, we consider a mode for which Xcntcr ^ 1- The second 
solution has the asymptotic form A(x) 2± ln(4/a;) for x <C 1. Using this and 
matching the solution at x = x^nter we get the properly matched mode, 
inside the Hubble radius, to be 



5{x) — Xguter 



1 + 2 



\ ^^enter / 



:i + Y)-2Xe\terln(-). (56) 
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During the radiation dominated phase — that is, till a ^ Oeq, x ^ I — this 
mode can grow by a factor 



d{x ~ 1) 
•^(^^enter) 



-^—5{x ~ 1) ^ 5In (^—) 

•^cntcr V S^enter / 

51n(^)=|ln(^). 
\ Center / 2 \ tenter / 



(57) 



Since the time tenter for a mode with wavelength A is fixed by the condition 

AOenter OC At^^^er - dn (tenter) OC tenter, it foUoWS that A OC tenter- HenCC, 

4nter V A ^ 3 V M ^ ^ ^ 

for a mode with wavelength A ^ Aeq. [Here, M is the mass contained in a 
sphere of radius (A/2); sec equation (9).] The growth in the radiation dom- 
inated phase, therefore, is logarithmic. Notice that the matching procedure 
has brought in an amplification factor which depends on the wavelength. 

In the discussion above, we have assumed that $7 = 1, which is a valid 
assumption in the early phases of the universe. However, during the later 
stages of evolution in a matter dominated phase, we have to take into 
account the actual value of n and solve equation (43). This can be done 
along the following lines. 

Let p(t) be a solution to the background Fricdmann model dominated 
by pressureless dust. Consider now the function pi(t) = p{t + T) where r is 
some constant. Since the Friedmann equations contain t only through the 
derivative, pi{t) is also a valid solution. If we now take r to be small, then 
[pi(t) — p(t)] will be a small perturbation to the density. The corresponding 
density contrast is 

^' Pit) Pit) dt ^ ^ 

where the last relation follows from the fact that p oc and H{t) = (d/a). 
Since r is a constant, it follows that H{t) is a solution to be the perturba- 
tion equation. [This curious fact, of course, can be verified directly: From 
the equations describing the Friedmann model, it follows that H + = 
(— 47rGp/3j. Differentiating this relation and using p = —3Hp we immedi- 
ately get H + 2HH —AirGpH = 0. Thus H satisfies the same equation as 

Since H = -H"^ - (47rGp/3), we know that < 0; that is, H \s a. 
decreasing function of time, and the solution J = = 5^ is a decaying 
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mode. The growing solution (S = 6g) can be again found by using the fact 
that, for any two hnearly independent solutions of the equation (42), the 
Wronskian {dgSj, —SdSg) has a value a~^. This implies that 

Hit) I -^g^. (60) 



a^j ^' J a-^m{t)' 

Thus we see that the H{t) of the background spacetime allows one to com- 
pletely determine the evolution of density contrast. 

It is more convenient to express this result in terms of the redshift z. 
For a universe with arbitrary fi, we have the relations 

a(z) = ao(l + z)-S H{z) = Ho{l + z){l + ^zf/'^ (61) 

and 

HQdt = -{l + z)-'^(l + nz)-^dz. (62) 
Taking 5d = H{z), we get 

Sg = 5d{z)ja-%\z){^f^dz 

roo 3 

= {aQHQ)-'^{l + z){l + ^zf''^j dx{l + x)-'^{l + ^x)-^. (63) 

This integral can be expressed in terms of elementary functions: 

_ l + 2^^ + 3^^z 3 n(l + z)(l + ^^z)V^ r (l + Qz)V^ + (l-r?)^/^ " 

^ ~ (1-Q)2 2 (1 - 0)5/2 ^ [(l + J^z)V2-(l-J^)l/2 

(64) 

Thus 6g{z) for an arbitrary can be given in closed form. The solution in 
(64) is not normalized in any manner; normalization can be achieved by 
multiplying 8g by some constant depending on the context. 

For large z (i.e., early times), 5g oc z^^. This is to be expected because 
for large z, the curvature term can be ignored and the Priedmann universe 
can be approximated as a = 1 model. [The large z expansion of the 
logarithm in (64) has to be taken upto 0{z~^/'^) to get the correct result; it 
is easier to obtain the asymptotic form directly from the integral in (63)]. 
For n <C 1, one can see that 5g ~ constant for z <C This is the 

curvature dominated phase, in which the growth of perturbations is halted 
by rapid expansion. 

We have thus obtained the complete evolutionary sequence for a per- 
turbation in the linear theory, which is shown in figure 1. This result can be 
conveniently summarized in terms of a quantity called 'transfer function' 
which we shall now describe. 
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Figure 1. Schematic figure showing the growth of linear perturbations in dark matter. 
The perturbation grows as before entering the Hubble radius when relativistic theory 
is required. During the radiation dominated phase it grows only as Ina and during the 
matter dominated phase it grows as a. In case the universe become dominated by curva- 
ture or background energy density, the perturbations do not grow significantly after that 
epoch. 

5. Transfer function 

If 5(t,x) <^ 1, then one can describe tlie evolution of x) by linear 
perturbation theory, in which each mode 5k (i) will evolve independently 
and we can write 

5k(t) = rk(i,ti)(5k(ii) (65) 

where T]^{t, ti) depends only on the dynamics and not on the initial condi- 
tions. We shall now determine the form of Tk{t, U). 

Let 5x{ti) denote the amplitude of the dark matter perturbation corre- 
sponding to some wavelength A at the initial instant t^. To each A, we can 
associate a wavenumber k oc A~^ and a mass M cx A^; accordingly, we may 
label the perturbation as SM{t) or 5k{t), as well, with the scalings M ~ A^, 
k ~ A~^. We are interested in the value of Sx{t) at some t ^ tdcc- 

To begin with, consider the modes which enter the Hubble radius in 
the radiation dominated phase; their growth is suppressed in the radiation 
dominated phase by the rapid expansion of the universe; therefore, they 
do not grow significantly until t = teq, giving 6x{teq) = (tenter) where 
L ~ 51n(Aeq/A) is a logarithmic factor determined in (58). After matter 
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begins to dominate, the amplitude of these modes grows in proportion to 
the scale factor a. Thus, 

5M{t) = L6M{tenter) ( — ) (for M < Meq). (66) 

Consider next the modes with Acq < A < Xh where A^^ = H^^{t) is the 
Hubble radius at the time t when we are studying the spectrum. These 
modes enter the Hubble radius in the matter dominated phase and grow 
proportional to a afterwards. So, 

Suit) = 6M{tenter). f^") (for Meq < M < Mh) (67) 
V ^^enter / 

which may be rewritten as 

Suit) = (^m(W) (^) (—] . (68) 

V O'enter J y ^eq J 

But notice that, since tenter is fixed by the condition Acenter o<^ Center 

Aienter> ^6 haVC ienter OC A^. Further (Oeq/aenter) = (teq/ienter)^*^^, giving 



\ '^ciitcr / 



Substituting (69) in (68), we get 



5M(t) = 5m(W) ('^y { — ) = <5m(W) { ^V^' { — ) . (70) 



A ; V"eqy V M . 

Comparing (70) and (66) we see that the mode which enters the Hubble 
radius after fcq has its amplitude decreased by a factor L~^M~^/^, com- 
pared to its original value. 

Finally, consider the modes with X > Xh which are still outside the 
Hubble radius at t and will enter the Hubble radius at some future time 
^cntcr > During the time interval (t, tenter), they will grow by a factor 

(aenter/a)- ThuS 

<5A(Wer) = m f^) (71) 



or 



^X{t) = <5a (tenter) ( = 5m(W) f^^^' ( — ] (A > Xh). 



VC^enter/ \ -^V^ 



(72) 
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[The last equality follows from the previous analysis]. Thus the behaviour 
of the modes is the same for the cases Aeq < \ < Xh and Xh < A; i.e. for 
all wavelengths A > Agq. Combining all these pieces of information, we can 
state the final result as follows: 

r /,x _/ ^^A (Center) (a/ Acq) (A < Agq) , . 

'^^W - \ 6x{tenter){a/a,^){X,JXf (Aeq < A) ^ 

or, equivalently 

X /,N /'^'^M(icntcr)(o/acq) (M < Mgq) 

' I 5M(tenter)(a/aeq)(Meq/M)2/3 (Meq < M). ^'^^ 

Thus the amplitude at late times is completely fixed by the amplitude of 
the modes when they enter the Hubble radius. 

In this approach, to determine S{x,t) or 6k{t) at time t, we need to 
know its exact space dependence (or k dependence) at some initial instant 
t = ti [eg. to determine S{t, x), we need to know 5{ti,x.)]. Often, we are not 
interested in the exact form of 6{t, x) but only in its "statistical properties" 
in the following sense: We may assume that, for sufficiently small U, each 
fourier mode dk{ti) was a Gaussian random variable with 

{SUti)6;{ti)) = (27r)3p(k, ti)5D{yi - P) (75) 

where P(k, tj) is the power spectrum of 6{ti,x) and < ••• > denotes an 
ensemble average. Then, 

{Sy,{t)6;{t)) = Tkit,u)T;{t,u){6Uu)6;{u)) 



(76) 

and the statistical nature of 6^ is preserved by evolution with the power 
spectrum evolving as 

P{k,t) = \T^{t,U)\''P{Kti). (77) 

It should be stressed that as far as linear evolution of perturbations are 
concerned the statistics of the perturbations is niaint aincd. For any random 
field one can define a power spectrum and study its evolution along the 
lines described below. In case of a gaussian random field with zero mean the 
power spectrum contains the complete information; in other cases the power 
spectrum will only provide partial information. This is the key difference 
between gaussian and other statistics. Some theories of structure formation 
describing the origin of initial perturbations predict the statistics of the 
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perturbations to be gaussian. Since this seems to be fairly natural we shall 
confine to this case in our discussion. 

A closely related quantity to the power spectrum is the two point cor- 
relation function, defined as 

6(x) = (J(x + y)6{y)) = J ^^(5k5;)e^k.(x+y)g-ip.y (^g) 

where < • • • > is the ensemble average. Using 

{5^5;) = {27TfPik)SD{k-p) (79) 

we get 

e5(x) = / ^P(k)e'''- (80) 

That is, the correlation function is the Fourier transform of the power 
spectrum. 

Our analysis can be used to determine the growth of P{k) or ^(x) as 
well. In practice, a more relevent quantity characterizing the density inho- 
mogeneity is A| = {k'^P{k) /2tt'^) where P{k) = |(5fcpis the power spectrum. 
Physically, A| represent the power in each logarithmic interval of k. From 
(73) we find that quantity behaves as 

A2 = (Center) (a/flcq)^ {f OT k^.^ < k) , . 

\Al{t,^tcr){a/a,^f{k/k,^)^ Uork<k,q). 

Let us next determine A| (tenter) if the initial power spectrum, when the 
mode was much larger than the Hubble radius, was a power law with A^ oc 
k^P{k) oc K^^^. This mode was growing as while it was bigger than 
the Hubble radius (in the radiation dominated phase). Hence A| (tenter) oc 
"^enter^"^^- the radiation dominated phase, we can relate aenter to A by 
noting that Aognter oc tenter OC agnterJ SO A OC Center OC Therefore, 

A2(tenter) OC a^nter^"+' OC (82) 

Using this in (81) we find that 

a2 _ / ^^(A;)/c"-i(a/aeq)2 (for A;eq < k) . . 

"^^-{k-^^ia/a^f (forfc<A;eq). ^^^^ 

This is the shape of the power spectrum for a > Ogq. It retains its initial 
primordial shape (A| oc fc"^+^) at very large scales (k < keq or A > Agq). 
At smaller scales, its amplitude is essentially reduced by four powers of k 
(from k'"'^^ to k'^~^). This arises because the small wavelength modes enter 
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the Hubble radius earlier on and their growth is suppressed more severely 

during the phase Center < O < Oeq- 

Note that the index n = 1 is special. In this case, (tenter) is inde- 
pendent of k and all the scales enter the Hubble radius with the same 
amplitude. The above analysis suggests that if n = 1, then all scales in 
the range keq < k will have nearly the same power except for the weak, 
logarithmic dependence through L'^{k). Small scales will have slightly more 
more power than the large scales due to this factor. 

There is another — completely different — reason because of which n = 
1 spectrum is special. If P{k) oc /c", the power spectrum for gravitational 
potential Pif{k) oc {P{k)/k'^) varies as P(p{k) oc fc""^. The power per loga- 
rithmic band in the gravitational potential varies as = (A;'^P(^(fe)/27r^) oc 
For n = 1, this is independnet of k and each logarithmic interval in k 
space contributes the same amount of power to the gravitational potential. 
Hence any fundamental physical process which is scale invariant will gen- 
erate a spectrum with n = 1. Thus observational verification of the index 
to n = 1 only verifies the fact that the fundamental process which led to 
the primordial fluctuations is scale invariant. 

Finally, we mention a few other related measures of inhomogeneity. 
Given a variable (5(x) we can smooth it over some scale by using win- 
dow functions W{x) of suitable radius and shape (We have suppressed the 
t dependence in the notation, writing 6{x,t) as 5(x)). Let the smoothed 
function be 

6w{^) = 1 5(x + y)W{y)(fy. (84) 
Fourier transforming 5vy(x), we find that 

If (5k is a Gaussian random variable, then Qk is also a Gaussian random 
variable. Clearly 5w{^) — which is obtained by adding several Gaussian 
random variables Qk — is also a Gaussian random variable. Therefore, to 
find the probability distribution of 5w (x) we only need to know the mean 
and variance of (5|y(x). These are, 

{Swi^))= J ^(<5k)W^^e^''-^ = 

((5^(x)) = J ^P(k)|Wk|^ ^ n'. (86) 
Hence the probability of 5w to have a value q at any location is given by 
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Note that this is independent of X, clS expected. 

A more interesting construct will be based on the following question: 
What is the probability that the value of 5\\/ at two points xi and X2 are 
qi and q2 ? Once we choose (xi,X2) the 6w {^i) , {^2) are correlated 
Gaussians with {6w (xi) Sw {'^2)) = (r) where r = xi — X2. The simul- 
taneous probability distribution for Sw{^i) = Qi and 5wi^2) = Q2 for two 
correlated Gaussians is given by: 



where 



'P[qi,q2] 



Q[qi,Q2] 



1 



27r/x2 



1/2 



exp-Q [51,^2] 



1-^2 



(88) 



(89) 



with A = [^Rif)/ fA- (This is easily verified by computing {qi),{q2) and 
{QiQj) explicitly). We can now ask: What is the probablility that both qi 
and q2 are high density peaks ? Such a question is particularly relevant since 
we may expect high density regions to be the locations of galaxy formation 
in the universe (see e.g. Kaiser, 1985). Then the correlation function of 
the galaxies will be the correlation between the high density peaks of the 
underlying gaussian random field. This is easily computed to be 



00 00 



P2 [qi > i'fJ,,q2 >iyij]= J dqi j dq2P[qi,q2] = Pi{q > i^fJ-) [1 + ^Ar)] 



UfJ, UfJ, 



(90) 

where ^i/(r) denotes the correlation function for regions with density which 
is I' times higher than the variance of the field. Explicit computation now 
gives 

00 00 

P2 oc y dti J dt2 exp{-^ ^372 (^1 + *i - 2^*1*2)} (91) 



U V 



This result can be expressed in terms of error function. An interesting 
special case in which this expression can be approximated occurs when 
^ ^ 1 and S> 1 though Av"^ is arbitrary. Then we get 



P2 = ^e"''' exp ^ {q > vjj.) exp (^z/^) 



(92) 



so that 



{f) = exp {j^^^^ ~ 1 = exp 



- 1 



(93) 
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In other words, the correlation function of high density peaks of a gaussian 
random field can be significantly higher than the correlation function of the 
underlying field. If we further assume that A S> 1 and Au^ ^ 1, 

then 



In this limit £,u{r) oc £,R{r) with the correlation increasing as i^^. 

A simple example of the window function arises in the following context. 
Consider the mass contained within a sphere of radius R centered at some 
point X in the universe. As we change x, keeping R constant, the mass 
enclosed by the sphere will vary randomly around a mean value Mq = 
{4:Tt/3)pbR^ where pB is the matter density of the background universe. 
The mean square fluctuation in this mass {{5M /M)\) is a good measure of 
the inhomogeneities present in the universe at the scale R. In this case, the 
window function is W{y) = 1 for |y| < R and zero otherwise. The variance 
in (86) becomes: 



This will be a useful statistic in many contexts. 

Another quantity which wc will use extensively in latter sections is the 
average value of the correlation function within a sphere of radius r, defined 




(94) 




(95) 



to be 




(96) 



Using 





and (96) we find that 




oo 




—P {k) [sin kr — kr cos kr] . 
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(98) 

A simple computation relates a^^^{R) to ^{x) and ^{x). We can show that 
and 




(100) 



Note that a^^^^ at R is determined entirely by ^(x) (or ^(x)) in the range 
< a; < 2R. (For a derivation, see Padmanabhan, 1996) 

The Gaussian nature of 6^ cannot be maintained if the evolution couples 
the modes for different values of k. Equation (36), which describes the 
evolution of 6k{t), shows that the modes do mix with each other as time 
goes on. Thus, in general, Gaussian nature of S^s cannot be maintained in 
the nonlinear epochs. 

6. Zeldovich approximation 

We shall next consider the evolution of perturbations in the nonlinear 
epochs. This is an intrinsically complex problem and the only exact proce- 
dure for studying it involves setting up large scale numerical simulations. 
Unfortunately numerical simulations tend to obscure the basic physics con- 
tained in the equations and essentially acts as a 'black box'. Hence it is 
worthwhile to analyse the nonlinear regime using some simple analytic ap- 
proximations in order to obtain insights into the problem. In sections 6 to 
8 and in section 11 we shall describe a series of such approximations with 
increasing degree of complexity. The first one — called Zeldovich approx- 
imation — is fairly simple and leads to an idea of the kind of structures 
which generically form in the universe. This approximation, however, is not 
of much use for more detailed work. The second and third approximations 
described in sections 7 and 8 arc more powerful and allow the modeling of 
the universe based on the evolution of the initially over dense region. Finally 
we discuss in section 11a fairly sophisticated approach involving nonlinear 
scaling relations which are present in the dynamics of gravitational clus- 
tering. In between the discussion of these approximations, we also describe 
some useful procedures which can be adopted to answer questions that are 
directly relevant to structure formation in sections 9 and 10. 

A useful insight into the nature of linear perturbation theory (as well as 
nonlinear clustering) can be obtained by examining the nature of particle 
trajectories which lead to the growth of the density contrast 5l(o) oc a. To 
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determine the particle trajectories corresponding to the hnear hmit, let us 
start by writing the trajectories in the form 

XT(a,q) = q + L(a,q) (101) 

where q is the Lagrangian coordinate (indicating the original postion of 
the particle) and L(a,q) is the displacement. The corresponding fourier 
transform of the density contrast is given by the general expression 

(5(a,k) = J (i3qg-ik.q-ik.L(a,q) _ (27r)3<5Dirac M (102) 

In the linear regime, we expect the particles to have moved very little and 
hence we can expand the integrand in the above equation in a Taylor series 
in (k • L). This gives, to the lowest order, 

d{a, ^) = - f d^fl e-'^ 'iiik ■L{a,q)) = - J (fq e-'^ '^ ( ■ L) (103) 

showing that 5{a, k) is Fourier transform of — Vq.L(a, q). This allows us to 
identify V • L(a, q) with the original density contrast in real space —6{a, q). 
Using the Poisson equation (for a $7 = 1, which is assumed for simplicity) 
we can write 6{a, q) as a divergence; that is 

V • L(a, q) = -6{a, q) = —H^^aV ■ (V</)) (104) 

which, in turn shows that a consistent set of displacements that will lead 
to (5(a) oc a is given by 

L(a, q) = -(VV')a = au(q); V = m)Ho^<t> (105) 

The trajectories in this limit are, therefore, linear in a: 

XT(a, q) = q + au(q) (106) 

An useful approximation to describe the quasilinear stages of clustering 
is obtained by using the trajectory in (106) as an ansatz valid even at 
quasilinear epochs. In this approximation, called Zeldovich approximation, 
the proper Eulerian position r of a particle is related to its Lagrangian 
position q by 

r{t) = a{t)yi{t) = a(t)[q + a(t)u(q)] (107) 

where x(t) is the comoving Eulerian coordinate. This relation in (106) gives 
the comoving position (x) and proper position (r) of a particle at time t, 
given that at some time in the past it had the comoving position q. If 
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the initial, unperturbed, density is p (which is independent of q), then the 
conservation of mass implies that the perturbed density will be 

p(r, t)d^Y = pcfq. (108) 

Therefore 



p{r,t) = p 



1 -1 



p/a^ Pbit) 



det {dxj I dqi ) det {Sij + a{t) {duj /dqi)) 

(109) 

where we have set Pb{t) = [p/a?{t)\. Since u(q) is a gradient of a scalar 
function, the Jacobian in the denominator of (109) is the determinant of a 
real symmetric matrix. This matrix can be diagonolized at every point q, 
to yield a set of eigenvalues and principal function of q. If the 

eigenvalues of (duj/dqi) are [— Ai(q), — A2(q), — A3(q)] then the perturbed 
density is given by 

= (1 - a(t)Ai(q))(l - r(t)A2(q))(l - a(i)A3(q)) 

where q can be expressed as a function of r by solving (107). This expression 
describes the effect of deformation of an infinitesimal, cubical, volume (with 
the faces of the cube determined by the eigenvectors corresponding to A^) 
and the consequent change in the density. A positive A denotes collapse and 
negative A signals expansion. 

In a ovcrdcnsc region, the density will become infinite if one of the terms 
in brackets in the denominator of (110) becomes zero. In the generic case, 
these eigenvalues will be different from each other; so that we can take, say, 
Ai > A2 > A3. At any particular value of q the density will diverge for the 
first time when (1 — a(t)Ai) = 0; at this instant the material contained in 
a cube in the q space gets compressed to a sheet in the r space, along the 
principal axis corresponding to Ai. Thus sheetlike structures, or 'pancakes', 
will be the first nonlinear structures to form when gravitational instability 
amplifies density perturbations. 

The trajectories in Zeldovich approximation, given by (106) can be used 
in (41) to provide a closed integral equation for In this case, 

XT(q,a) = q + aVV^; xt = (^^^ Vip; = ^p'P (m) 

and, to the same order of accuracy, i?k in (38) becomes: 

J d^q (k ■ xt)^ e-'^<'^+^^ ^ J d3q(k ■ XT)^e-^'*-'i (112) 
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Substituting these expressions in (41) we find that the gravitational poten- 
tial is described by the closed integral equation: 

(113) 

This equation provides a powerful method for analysing non linear cluster- 
ing since estimating [A\^ — Bi^) by Zcldovich approximation has a very large 
domain of applicability (Padmanabhan, 1998). 

It is also possible to determine the power spectrum corresponding to 
these trajectories using our general formula 

P(k, a) = \S{k, a)|2 = J d^qd^q'e-^'^-^i-i') <Je-*-[^("''i)-^(''''i')]) (114) 

The ensemble averaging can be performed using the general result for gaus- 
sian random fields: 

(^e^^-^^ = cxp (^-kikja'^{V)/2^ (115) 

where cr*-' is the covariance matrix for the components V"' of a gaussian 
random field. This quantity can be expressed in terms of the power spec- 
trum Piik) in the linear theory and a straightforward analysis gives (see, 
for e.g., Taylor and Hamilton, 1996) 

POO /■+! p -, 

P{k,a)=J i-Kq^dqj d^le'^'i^' eycp-k^ [F{q) + ti^qF'{q)\ (116) 

where 

The integrals, unfortunately, needs to be evaluated numerically except in 
the case of n = —2. In this case, we get 



A2(A;,a) = 



k^P 16 



a^k 



27r2 TT [1 + (2a2fc)2]2 



1 + 



37r 



a^k 



4 [1 + (2a2fc)2]i/2 



(118) 



which shows that A2 (x a? for small a but decays as a~2 at late times due to 
the dispersion of particles. Clearly, Zeldovich approximation breaks down 
beyond a particular epoch and is of limited validity. 
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7. Spherical approximation 

In the nonhnear regime — when (5^1 — it is not possible to solve equa- 
tion (36) exactly. Some progress, however, can be made if we assume that 
the trajectories are homogeneous; i.e. x(t, q) = /(t)q where f{t) is to be 
determined. In this case, the density contrast is 

= (2^)3,5D(k)[/-3-l] = (27r)35i)(k)5(t) (119) 

where we have defined 5{t) = [f^^'{t) — l] as the amplitude of the density 
contrast for the k = mode. It is now straightforward to compute A and 
B in (36). We have 

A = 4TrGpbS'^{t)[{2TrfSD0<i)] (120) 



and 

B 



I d^qik^afpe-'f^'^^i"^ = -/2^[(27r)3fe(/k)] 
^ [(27r)=^^z.(k)] (121) 



3(1 + 5) 

so that the equation (36) becomes 

S + 2^S = 47rGp6(l + S)S + (122) 

(This particular approach to spherical collapse model, which does not re- 
quire fluid equations is due to Padmanabhan 1998.) To understand what 
this equation means, let us consider, at some initial epoch tj, a spherical 
region of the universe which has a slight constant overdensity compared to 
the background. As the universe expands, the over dense region will expand 
more slowly compared to the background, will reach a maximum radius, 
contract and virialize to form a bound nonlinear system. Such a model is 
called "spherical top-hat" . For this spherical region of radius R{t) contain- 
ing dustlike matter of mass M in addition to other forms of energy densities, 
the density contrast for dust will be given by: 



l + S 



p 3M 1 2GM 



ait) 



R{t) 



[Note that, with this definition / oc (R/a).] Using this in (122) we can 
to obtain an equation for R{t) from the equation for 5; straight forward 
analysis gives 

GM AttG 

'r^ 
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This equation could have been written down "by inspection" using the 
relations 

R = -V</)tot; (Ptot = (pFRW + Scl> = -{d/2a)R^ - G5M/R. (125) 

Note that this equation is valid for perturbed "dust-like" matter in any 
background spacetime with density prest aiid pressure Prest contributed by 
the rest of the matter. Our homogeneous trajectories x(q, t) = /(t)q actu- 
ally describe the spherical top hat model. 

This model is particularly simple for the = 1, matter dominated 
universe, in which prest = Prest = and we have to solve the equation 

<^R GM 

This can be done by standard techniques and the final results for the evo- 
lution of a spherical overdense region can be summarized by the following 
relations: 

= 2^^^ " "^'^^ = 10^^^ " ^^^^^ 



,,9(e-sin^)2 ^^^^^ 
''W = ''>W 2(l-cos^)3 - 

The density can be expressed in terms of the redshift by using the relation 
{t/tif/^ = (1 + Zi){l + z)-^. This gives 

^^^^^{'i) (0-sin^)V3 = (3) (3) (0-sin Of/^' 



u yj 0111 w 



9 (^ - sin ef 
2(1- cos ef 

Given an initial density contrast 5i at redshift z^, these equations define 
(impHcitly) the function 5{z) for z > Zi. Equation (130) defines 9 in terms 
of z (implicitly); equation (131) gives the density contrast at that 0[z). 

For comparison, note that linear evolution gives the density contrast 5^ 
where 



Ph b l + z 5 ^ 



sin^)2/3. (132) 



We can estimate the accuracy of the linear theory by comparing 5{z) and 
5l{z). To begin with, for z » 1, we have <C 1 and we get 5{z) ~ 
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6l{z). When 6 = (7r/2), 5l = (3/5)(3/4)2/3 (7r/2 - l)^/^ = 0.341 while 
6 = (9/2)(7r/2 — 1)2 — 1 = 0.466; thus the actual density contrast is about 
40 percent higher. When 6 = {2tt/?,),5l = 0.568 and 5 = 1.01 ~ 1. If we 
interpret 5 = 1 as the transition point to nonlinearity, then such a transi- 
tion occurs at ^ = (27r/3), 8l — 0.57. Prom (130), we see that this occurs 
at the redshift (1 + zj) = 1.065^(1 + Zi) = (5o/0.57). 

The spherical region reaches the maximum radius of expansion at = vr. 
From our equations, we find that the redshift Zm, the proper radius of the 
shell rm and the average density contrast 5^ at 'turn-around' are: 



5 (5( 







3(37r/4)2/3 1.062' 

'^x f-p\ , ^ 97r2 
r-m. = — = 1 -l-d^ = — ~ 5.6. 

Soo \PhJm 16 



(133) 



The first equation gives the redshift at turn-around for a region, parametrized 
by the (hypothetical) linear density contrast 6q extrapolated to the present 
epoch. If, for example, 5i ~ 10"^ at zi ~ 10^, such a perturbation would 
have turned around at {1+Zm) — 5.7 or when Zm — 4.7. The second equation 
gives the maximum radius reached by the perturbation. The third equation 
shows that the region under consideration is nearly six times denser than 
the background universe, at turn-around. This corresponds to a density 
contrast of 5m ~ 4.6 which is definitely in the nonlinear regime. The linear 
evolution gives 5l = 1.063 at ^ = tt. 

After the spherical overdense region turns around it will continue to con- 
tract. Equation (129) suggests that at = 27r all the mass will collapse to 
a point. However, long before this happens, the approximation that matter 
is distributed in spherical shells and that random velocities of the par- 
ticles are small, (implicit in the assumption of homogeneous trajectories 
X = /(i)q) will break down. The collisionless (dark matter) component 
will relax to a configuration with radius Tvir, velocity dispersion v and den- 
sity PcoU- After virialization of the collapsed shell, the potential energy U 
and the kinetic energy K will be related by \U\ = 2K so that the total 
energy £ = U + K = ~K . At t = tm sdl the energy was in the form of 
potential energy. For a spherically symmetric system with constant density, 
£ ~ —3GM'^/5rrn- The 'virial velocity' v and the 'virial radius' ryir for the 
collapsing mass can be estimated by the equations: 

K^^ = -£ =^—; \U\ = ^-- = 2K = Mv\ (134) 

z <J' m 'Ji vir 
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We get: 

V = (6GM/5r„)^/2. ^^.^ ^ ^^/2. (135) 

The time taken for the fluctuation to reach virial equilibrium, tcow, is es- 
sentially the time corresponding to 6 = 27r. Prom equation (130), we find 
that the redshift at collapse, Zcou, is 

(l+^coii) = (2^^^3(^4)2/3 = 0.365.(1+..) = 0.63(1+.^) = (136) 

The density of the collapsed object can also be determined fairly easily. 
Since rvir = (?'m/2), the mean density of the collapsed object is pcoii = ^Pm 
where pm is the density of the object at turn-around. 

We have, pm = 5.Qpb{tm) and pbitm) = (1 + ^mf (1 + -2:coll)"^P6(icoll)- 
Combining these relations, we get 

PcoU - 'i^Pm ^ ^.8pb{tm) ^ 170p6(tcoll) - 170po(l + ^cou)^ (137) 

where po is the present cosmological density. This result determines Pcoii in 
terms of the redshift of formation of a bound object. Once the system has 
virialized, its density and size does not change. Since pb oc a~^, the density 
contrast 6 increases as for t > tcoii- 

This approach can be easily generalised to describe the situation in 
which the initial density profile is given by p{r.i). Given an initial density 
profile Pi(r), we can calculate the mass M{ri) and energy E{ri) of each 
shell labelled by the initial radius r.. In spherically symmetric evolution, 
M and E are conserved and each shell will be described by equation (126). 
Assuming that the average density contrast 5i{ri) decreases with r., the 
shells will never cross during the evolution. Each shell will evolve in accor- 
dance with the equations (127), (128) with 5i replaced by the mean initial 
density contrast 5j(rj) characterising the shell of initial radius r^. Equation 
(129) gives the mean density inside each of the shells from which the density 
profile can be computed at any given instant. 

A simple example for this case corresponds to a scale invariant situation 
in which E{M) is a power law. If the energy of a shell containing mass M 
is taken to be 

E{M) = Eo l^—J < 0, (138) 



then the turn-around radius and turn-around time are given by 

GM _ GMo ( M 

E{M) ~ E^ yWo 



(M)^-^^-^^fily^' (139) 
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To avoid shell crossing, we must have e > so that outer shells with more 
mass turn around at later times. In such a scenario, the inner shells expand, 
turn around, collapse and virialize first and the virialization proceeds pro- 
gressively to outer shells. We shall assume that each virialized shell settles 
down to a final radius which is a fixed fraction of the maximum radius. 
Then the density in the virialized part will scale as (M/r^) where M is the 
mass contained inside a shell whose turn-around radius is r. Using (139) to 
relate the turn-around radius and mass, we find that 

p^r) oc ''^ oc r3/(i+3^)r-3 oc r-'^/(^+^^\ (141) 

Two special cases of this scaling relation are worth mentioning: (i) If the 
energy of each shell is dominated by a central mass m located at the origin, 
then E oc Gm/r oc M~^/^. In that case, e = 1 and the density profile of 
virialized region falls as r~^/^. The situation corresponds to a accretion on 
to a massive object (ii) If e = 2/3 then the binding energy E is the same for 
all shells. Then we get p oc which corresponds to an isothermal sphere. 

The spherical model can be easily generalised for the set of trajectories 
with x"(t, q) = f"'^{t)qi) (Padmanabhan 1998) In this case, it is convenient 
to decompose the derivative of the velocity daU\, = fab into shear aab, 
rotation O'^ and expansion 9 by writing 

fab = (Tab + eabc^" + ^<5a&^. (142) 

where a^b is the symmetric traceless part of /ab; the e^.^c^'^ is the antisym- 
metric part and {l/3)5ah& is the trace. In this case, (122) gets generalised 
to: 

4 

6 + 2-6 = 47rGp6(l + 6)6 + + 0^(1 + 6)ia^ - 2^"^) (143) 

a 6 [1 + 0) 

where = (J°'^CFab ^-nd = . From the last term on the right 

hand side we see that shear contributes positively to 6 while rotation 
contributes negatively. Thus shear helps growth of inhomogenitics while ro- 
tation works against it. To see this explicitly, we again introduce a function 
R{t) by the definition 

^ , 9GMt2 a3 ^^^^^ 
l + ^=^-M^ (144) 
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where M and n are constants. Using this relation between 6 and R{t), 
equation (143) can be converted into the following equation for R{t) 

R = -^-^a'^a'-2n')R (145) 

where the first term represents the gravitational attraction diic to the mass 
inside a sphere of radius R and the second gives the effect of the shear and 
angular momentum. We shall now see how an improved spherical collapse 
model can be constructed with this term. 

8. Improved spherical collapse model 

In the spherical collapse model (SCM, for short) each spherical shell ex- 
pands at a progressively slower rate against the self-gravity of the system, 
reaches a maximum radius and then collapses under its own gravity, with a 
steadily increasing density contrast. The maximum radius, Rmax = Ri/^i, 
achieved by the shell, occurs at a density contrast 5 = (Qvr^/IG) — 1 ~ 4.6, 
which is in the "quasi-linear" regime. In the case of a perfectly spherical 
system, there exists no mechanism to halt the infall, which proceeds in- 
exorably towards a singularity, with all the mass of the system collapsing 
to a single point. Thus, the fate of the shell is to collapse to zero radius 
at ^ = 27r with an infinite density contrast; this is, of course, physically 
unacceptable. 

In real systems, however, the implicit assumptions that (i) matter is 
distributed in spherical shells and (ii) the non-radial components of the 
velocities of the particles are small, will break down long before infinite 
densities are reached. Instead, we expect the collisionless dark matter to 
reach virial equilibrium. After virialization, \U\ = 2K, where U and K 
are, respectively, the potential and kinetic energies; the virial radius can be 
easily computed to be half the maximum radius reached by the system. 

The virialization argument is clearly physically well-motivated for real 
systems. However, as mentioned earlier, there exists no mechanism in the 
standard SCM to bring about this virialization; hence, one has to introduce 
by hand the assumption that, as the shell collapses and reaches a particular 
radius, say Rmaxf^j the collapse is halted and the shell remains at this 
radius thereafter. This arbitrary introduction of virialization is clearly one 
of the major drawbacks of the standard SCM and takes away its predictive 
power in the later stages of evolution. Wc shall now see how the retention 
of the angular momentum term in equation (145) can serve to stabilize the 
collapse of the system, thereby allowing us to model the evolution towards 
r^ir = Rmaxf^ smoothly. (Engineer etal, 1998) 

At this point, it is important to note a somewhat subtle aspect of our 
generalisation. The original equations are clearly Eulerian in nature: i.e. 
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the time derivatives give the temporal variation of the quantities at a fixed 
point in space. However, the time derivatives in equation (143), for the 
density contrast 5, arc of a different kind. Here, the observer is moving 
with the fluid element and hence, in this, Lagrangian case, the variation 
in density contrast seen by the observer has, along with the intrinsic time 
variation, a component which arises as a consequence of his being at dif- 
ferent locations in space at different instants of time. When the 5 equation 
is converted into an equation for the function R{t), the Lagrangian pic- 
ture is retained; in SCM, we can interpret R{t) as the radius of a spherical 
shell, co-moving with the observer. The mass M within each shell remains 
constant in the absence of shell crossing and the entire formalism is well 
defined. The physical identification of R is, however, not so clear in the case 
where the shear and rotation terms are retained, as these terms break the 
spherical symmetry of the system. We will nevertheless continue to think 
of R as the "efl'ective shell radius" in this situation, defined by equation 
(144) governing its evolution. Of course, there is no such ambiguity in the 
mathematical definition of R in this formalism. This is equivalent to taking 
R^ as proportional to the volume of a region defined by the location of a 
set of mass points. 

We now return to equation (143), and recast the equation into a form 
more suitable for analysis. Using logarithmic variables, D^c = In (1 + 5) 
and a = Ina, equation (143) can be written in the form (the subscript 'SC 
stands for 'Spherical Collapse') 

d^Dsc 1 fdPscy j IdDsc 
do?' ?)\ da ) 2 da 

^ [exp{Dsc) - 1] + a'(a2 - 2n^) (146) 

where a takes the role of time coordinate. It is also convenient to introduce 
the quantity, S, defined by 

S = a^{a^ - 2^2) (147) 

which we shall hereafter call the "virialization term" . The consequences of 
the retention of the virialization term are easy to describe qualitatively. 
We expect the evolution of an initially spherical shell to proceed along the 
lines of the standard SCM in the initial stages, when any deviations from 
spherical symmetry, present in the initial conditions, are small. However, 
once the maximum radius is reached and the shell recollapses, these small 
deviations are amplified by a positive feedback mechanism. To understand 
this, we note that all particles in a given spherical shell are equivalent due to 
the spherical symmetry of the system. This implies that the motion of any 
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particle, in a specific shell, can be considered representative of the motion 
of the shell as a whole. Hence, the behaviour of the shell radius can be 
understood by an analysis of the motion of a single particle. The equation 
of motion of a particle in an expanding universe can be written as 

X, + 2^X. = -^ (148) 

where a{t) is the expansion factor of the locally overdense "universe" . The 
Xj term acts as a damping force when it is positive; i.e. while the back- 
ground is expanding. However, when the overdense region reaches the point 
of maximum expansion and turns around, this term becomes negative, act- 
ing like a negative damping term, thereby amplifying any deviations from 
spherical symmetry which might have been initially present. Non-radial 
components of velocities build up, leading to a randomization of velocities 
which finally results in a virialised structure, with the mean relative veloc- 
ity between any two particles balanced by the Hubble flow. It must be kept 
in mind, however, that the introduction of the virialization term changes 
the behaviour of the solution in a global sense and it is not strictly correct 
to say that this term starts to play a role only after recollapse, with the 
evolution proceeding along the lines of the standard SCM until then. It 
is nevertheless reasonable to expect that, at early times when the term is 
small, the system will evolve as standard SCM to reach a maximum radius, 
but will fall back smoothly to a constant size later on. 

Equation (143) is actually valid for any fluid system and the virialization 
term, S, is, in general, a function of a and x, since the derivatives in equation 
(143) are total time derivatives, which, for an expanding Universe, contain 
partial derivatives with respect to both x and t separately. Even in the case 
of displacements with = f°'^{t)qh, the one equation (143) cannot uniquely 
determine all the components of f°'^{t). Handling this equation exactly will 
take us back to the full non-linear equations and, of course, no progress 
can be made. Instead, we will make the ansatz that the virialization term 
depends on t and x only through 5(t,x): 

S(a,x) = S((5(a,x)) = 5(L>sc) (149) 

In other words, S" is a function of the density contrast alone. This ansatz 
seems well motivated because the density contrast, 5, can be used to char- 
acterize the SCM at any point in its evolution and one might expect the 
virialization term to be a function only of the system's state, at least to the 
lowest order. Further, the results obtained with this assumption appear to 
be sensible and may be treated as a test of the ansatz in its own framework. 
To proceed further systematically, we define a function hsc by the relation 
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^^3,.e (150, 

For consistency, we shall assume the ansatz hsc{ci,x.) = hsc [S{a,x)]. The 
definition of hsc allows us to write equation (146) as 

= hsc-^ + 2 i^MDsc) - 1] + 3 (151) 

Dividing (151) by (150), we obtain the following equation for the function 
^sc(-Dsc) 

If we know the form of either hsc{Dsc) or S{Dsc), this equation allows us 
to determine the other. Then, using equation (150), one can determine -Dsc- 
Thus, our modification of the standard SCM essentially involves providing 
the form of 5*30 (-Dsc) or /isc(-Dsc)- We shall now discuss several features 
of such a modelling in order to arrive at a suitable form. 

The behaviour of hsciDsc) can be qualitatively understood from our 
knowledge of the behaviour of S with time. In the linear regime {6 <^ 1), 
we know that S grows linearly with a; hence ^sc increases with Dsc- At 
the extreme non-linear end (5^1), the system "virializes" , i.e. the proper 
radius and the density of the system become constant. On the other hand, 
the density p^, of the background, falls like (or a~^) in a fiat, dust- 
dominated universe. The density contrast is defined by (5 = {p/pb—l) ^ p/pb 
(for (5 S> 1) and hence 

(5 oc oc (153) 

in the non-linear limit. Equation (150) then implies that hsci^) tends to 
unity for 5 3> 1. Thus, we expect that /isc(-C^sc) will start with a value 
far less than unity, grow, reach a maximum a little greater than one and 
then smoothly fall back to unity. [A more general situation discussed in the 
literature corresponds to h ^ constant as S oo, though the asymptotic 
value of h is not necessarily unity. Our discussion can be generalised to this 
case.] 

This behaviour of the hsc function can be given another useful inter- 
pretation whenever the density contrast has a monotonically decreasing 
relationship with the scale, x, with small x implying large 6 and vice- versa. 
Then, if we use a local power law approximation d cx x~'^ for 6^1 with 
some n > 0, we have Dsc ln(x^^) and 

dDsc dln{l) xa v 

hsc ^ — ; — — — oc — oc (154) 

da a In a ax ax 
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where v = ax denotes the mean relative velocity. Thus, hsc is proportional 
to the ratio of the relative peculiar velocity to the Hubble velocity. We 
know that this ratio is small in the linear regime (where the Hubble flow 
is dominant) and later increases, reaches a maximum and finally falls back 
to unity with the formation of a stable structure; this is another argument 
leading to the same qualitative behaviour of the hsc function. 

Note that, in standard SCM (for which S = 0), equation (152) reduces to 

3hscP^ = hlc-^ + - (155) 

The presence of the linear term in 6 on the RHS of the above equation causes 
hsc to increase with 6, with hsc 5^^"^ for (5 3> 1. If virialization is imposed 
as an ad hoc condition, then hsc should fall back to unity discontinuously 
— which is clearly unphysical; the form of S{6) must hence be chosen so as 
to ensure a smooth transition in /isc(<^) from one regime to another. [As an 
aside, we remark that S{5) can be reinterpreted to include the lowest order 
contributions arising from shell crossing, multi-streaming, etc., besides the 
shear and angular momentum terms, i.e. it contains all effects leading to 
virialization of the system; see S. Engineer, etal, 1998] 

We will now derive an approximate functional form for the virialization 
function from physically well-motivated arguments. If the virialization term 
is retained in equation (145), we have 

d^R GM H^R 

where H = d/a. Let us first consider the late time behaviour of the system. 
When virialization occurs, it seems reasonable to assume that R — > constant 
and R^ 0. This implies that, for large density contrasts, 

^--^ (157) 

Using H = d/a = {2/3t), and equation (144) 

_ 2>7GAlt^ 3 3 ^ , ^ . , . _ 

5~ 4R^ = ~2^^'^^^'^~2^ ('^»1) (15^) 

Thus, the "virialization" term tends to a value of (—35/2) in the non-linear 
regime, when stable structures have formed. This asymptotic form for S{6) 
is, however, insufficient to model its behaviour over the larger range of 
density contrast (especially the quasi-linear regime) which is of interest 
to us. Since S{S) tends to the above asymptotic form at late times, the 
residual part, i.e. the part that remains after the asymptotic value has 
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been subtracted away, can be expanded in a Taylor series in (1/(5) without 
any loss of generality. Retaining the first two terms of expansion, we write 
the complete virialization term as 

S{6) = -l{l + 6)-j + ^ + 0{6-') (159) 
Replacing for S{d) in equation (146), we obtain, for 5 S> 1 

3M^-/.sc + ^ + 2=-y + ^ (160) 

[It can be easily demonstrated that the first order term in the Taylor series 
is alone insufficient to model the turnaround behaviour of the h function. 
We will hence include the next higher order term and use the form in 
equation (159) for the virialization term. The signs are chosen for future 
convenience, since it will turn out that both A and B are greater than 
zero.] In fact, for sufficiently large 6, the evolution depends only on the 
combination q = (B/A^). Equation (156) can be now written as 



" GM 4i? 
R = 



B? 27^2 



27GMt^ A B 
+ 



(161) 



Using 5 = 9GMt^ /2R"' and B = qA^ we may express equation (161) com- 
pletely in terms of R and t. We now rescale R and t in the form R = rviry{x) 
and t = f3x, where r^ir is the final virialised radius [i.e. R — r^r for t — 00], 
and 0^ = (8/3^)(^/GM)r^j^, to obtain the following equation for y{x) 

We can integrate this equation to find a form for yq(.'r) (where Dqix) is 
the function y(.T) for a specific value of q) using the physically motivated 
boundary conditions y = \ and y' = as a; ^ 00, which is simply an 
expression of the fact that the system reaches the virial radius r^r and 
remains at that radius asymptotically. The results of numerical integration 
of this equation for a range of q values are shown in figure (2) . As expected 
on physical grounds, the function has a maximum and gracefully decreases 
to unity for large values of x [the behaviour of y{x) near x = is irrelevant 
since the original equation is valid only for b > 1, at least]. For a given 
value of g, it is possible to find the value Xc at which the function reaches 
its maximum, as well as the ratio ymax = Rmax/fvir- The time, tmaxi at 
which the system will reach the maximum radius is related to Xc by the 
relation tmax = I^Xc = io(l + Zmax)~^^'^i where = 2/{3Hq) is the present 
age of the universe and Zj^ax 

is the redshift at which the system turns 
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Figure 2. The figure shows the function yq{x) for some values of q. The x axis has scaled 
time, X and the y axis is the scaled radius y. 

around. Fig ure (3) shows the variation of and Umax — {-^max / fvir^ 
different values of q. The entire evolution of the system in the modified 
spherical collapse model (MSCM) can be expressed in terms of 



where /? = ito/xc){l + Zmax) ^Z^- 

In SCM, the conventional value used for [r-uir / Rmax) is (1/2), which is 
obtained by enforcing the virial condition that \U\ = 2K, where U is the 
gravitational potential energy and K is the kinetic energy. It must be kept 
in mind, however, that the ratio {vmr / Rmax) is not really constrained to be 
precisely (1/2) since the actual value will depend on the final density profile 
and the precise definitions used for these radii. While we expect it to be 
around 0.5, some amount of variation, say between 0.25 and 0.75, cannot 
be ruled out theoretically. 

Figure (3) shows the parameter {Rmax/'Tvir)-, plotted as a function of 
q = B/A^ (dashed line), obtained by numerical integration of equation 
(156) with the ansatz (159). The solid line gives the dependence of Xc (or 
equivalently tmax) on the value of q. It can be seen that one can obtain a 
suitable value for the {vvir/Rmax) ratio by choosing a suitable value for q 
and vice versa. 



R{t) = r^ir yq{t//3) 



(163) 
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Figure 3. The figure shows the parameters [Rmax/rvir) (broken line) and Xc (solid 
line) as a function of q = B/A?. This clearly demonstrates that the single parameter 
description of the virialization term is constrained by the value that is chosen for the 

ratio ?"i?zr /^max ■ 



Using equation (150) and the definition 6 cx t /R , we obtain 

hsc{x) = l-^-^ (164) 
2y ax 

which gives the form of h^cix) for a given value of q; this, in turn, deter- 
mines the function yq{x). Since 5 can be expressed in terms of x, y and 
Xc as (5 = (97r^/2a:^)x^/y^, this ahows us to implicitly obtain a form for 
hsc{^)i determined only by the value of q. 

It is possible to determine the best-fit value for q by comparing these 
results with simulations. This is best done by comparing the form of hsc{x). 
Such an analysis gives q = 0.02. (see S. Engineer, etal., 1998) Figure 
(4) shows the plot of scaled radius yq{x) vs x, obtained by integrating 
equation(162), with q = 0.02. The figure also shows an accurate fit (dashed 
line) to this solution of the form 

with a = —3.6, 6 = 53 and c = —12. This fit, along with values for r^r 
and Zmax-, completely specifies our model through equation (163). It can be 
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Figure 4- The figure shows a plot of the scaled radius of the shell as a function of 
scaled time x (solid line) and the fitting formula yq = {x + ax^ + bx^)/{l + cx^ + bx^), 
with a = —3.6, & = 53 and c = — 12 (dashed line) (See text for discussion) 



observed that {rvir / Rmax) is approximately 0.65. It is interesting to note 
that the value obtained for the {vyir / Rmax) ratio is not very widely off the 
usual value of 0.5 used in the standard spherical collapse model, in spite of 
the fact that no constraint was imposed on this value, ah initio, in arriving 
at this result. Finally, figure (5) compares the non-linear density contrast 
in the modified SCM (dashed line) with that in the standard SCM (solid 
line), by plotting both against the linearly extrapolated density contrast, 
5l- It can be seen (for a given system with the same Zmax and Vmr) that, 
at the epoch where the standard SCM model has a singular behaviour 
{5l ~ 1.686), our model has a smooth behaviour with 5 ^ 110 (the value is 
not very sensitive to the exact value of q). This is not widely off from the 
value usually obtained from the ad hoc procedure applied in the standard 
spherical collapse model. In a way, this explains the unreasonable effective- 
ness of standard SCM in the study of non- linear clustering. Figure (5) also 
shows a comparison between the standard SCM and the MSCM in terms 
of 6 values in the MSCM at two important epochs, indicated by vertical 
arrows, (i) When R = Rmax/^ in the SCM, i.e. the epoch at which the 
SCM virializes, (5(MSCM) ~ 83. (n) When the SCM hits the singularity, 
{Sl ~ 1.6865), 5(MSCM) - 110. 
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Figure 5. The figure shows the non-linear density contrast in the SCM (sohd line) 
and in the modified SCM (dashed line) , plotted against the linearly extrapolated density 
contrast Sl (discussion in text). 



9. Mass functions and abundances 

The description developed so far can also be used to address an important 
question: What fraction of the matter in the universe has formed bound 
structures at any given epoch and what is the distribution in mass of these 
bound structures? We shall now describe a simple approach which answers 
these questions. (Press and Schechter,1974) 

Gravitationally bound objects in the universe, like galaxies, span a large 
dynamic range in mass. Let f{M)dM be the number density of bound ob- 
jects in the mass range (M, M + dM) [usually called the "mass function"] 
and let F(M) be the number density of objects with masses greater than 
M. Since the formation of gravitationally bound objects is an inherently 
nonlinear process, it might seem that the linear theory cannot be used to 
determine F{M). This, however, is not entirely true. In any one realization 
of the linear density field 5i?(x), filtered using a window function of scale 
R, there will be regions with high density [i.e. regions with 5r > 5c where 
Sc is some critical value slightly greater than unity, say]. It seems reason- 
able to assume that such regions will eventually condense out as bound 
objects. Though the dynamics of that region will be nonlinear, the process 
of condensation is unlikely to change the mass contained in that region sig- 
nificantly. Therefore, if we can estimate the mean number of regions with 
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> (5c in a Gaussian random field, we will be able to determine F{M). 

One way of achieving this is as follows: Let us consider a density field 
(5/j(x) smoothed by a window function W/j of scale radius R. As a first 
approximation, we may assume that the region with 5{R,t) > 5c (when 
smoothed on the scale R at time U) will form a gravitationally bound object 
with mass M oc pR^ by the time t. The precise form of the M — R relation 
depends on the window function used; for a step function M = (47r/3) 
'pR^, while for a Gaussian M = (27r)^^'^'pR^. Here dc is a critical value 
for the density contrast which has to be supplied by theory. For example, 
dc — 1.68 in spherical collapse model. Since 5 oc 

^2/3 for a n = 1 universe, 
the probability for the region to form a bound structure at t is the same as 
the probability 6 > Sc{ti/t)'^^^ at some early epoch ti. This probability can 
be easily estimated since at sufficiently early ti, the system is described by 
a gaussian random field. Hence fraction of bound objects with mass greater 
than M will be 



F{M) = r P{S,R,ti)dS = ^-^ 

JSc{t,ti) V27T<7{R,t 



1 . / Sc{t,t, 

-erfc 




2 \V2a{R,ti 

where erfc(x) is the complementary error function. The mass function f{M) 
is just (dF/dM); the (comoving) number density N[M, t) can be found by 
dividing this expression by (M/'p). Carrying out these operations we get 



N(M,tm = - [fi) [^) [i) [-^) exp ^-^ j m. (167) 

Given the power spectrum \6k\'^ and a window function Wr one can explic- 
itly compute the right hand side of this expression. 

There is, however, one fundamental difficulty with the equation (166). 
The integral of /(M) over all M should give unity; but it is easy to see 
that, for the expression in (166), 

/ f{M)dM = dF = -. (168) 
Jo Jo ^ 

This arises because we have not taken into account the underdense regions 
correctly. To see the origin of this difficulty more clearly, consider the in- 
terpretation of (166). If a point in space has S > Sc when filtered at scale 
R, then that point should correspond to a system with mass greater than 
M{R); this is taken care of correctly by equation (166). However, consider 
those points which have 6 < Sc under this filtering. There is a non-zero 
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probability that such a point wiU have 6 > 6c when the density field is 
filtered with a radius Ri > R. Therefore, to be consistent with the inter- 
pretation in (166), such points should also correspond to a region with mass 
greater than M. But (166) ignores these points completely and thus under- 
estimates F{M) [by a factor (1/2)]. To correct this, we shall 'renormalise' 
the result by multiplying it by a factor 2. Then 

or 

The quantity a here refers to the linearly extrapolated density ai', the 
subscript L is omitted to simplify notation. The corresponding result for 
F{M) is also larger by factor two: 

(171) 

where o"o(M) is the linearly extrapolated density contrast today and we 
have used the fact (Tl(M, z) oc {1 + z)~^. Note that, by definition, F{M, z) 
gives the Jl contributed by the collapsed objects with mass larger than M 
at redshift z; equation (171) shows that this can be calculated given only 
the linearly extrapolated (Jo(M). The top panel of figure (6) gives 0,(M) as 
a function of ao{M) for different z, and the observed abundance of different 
structures in the universe. The six different curves from top to bottom are 
for z = 0,1, 2, 3, 4, 5. The dashed line on the z = curve gives the observed 
abundance of clusters; the trapezoidal region between z = 2 and z = 3 
is based on abundance of damped lyman alpha systems; the line between 
z = 2 and 2; = 4 is a lower bound on quasar abundance. 

As an example, let us consider the abundance of Abell clusters. Let 
the mass of Abell clusters to be M = 5 x lO^^aM© where a quantifies 
our uncertainty in the observtion. Similarly, we take the abundance to be 
^ = 4 X 10~^/3/i^Mpc~^ with (3 quantifying the corresponding uncertainty. 
The contribution of the Abell clusters to the density of the universe is 

F = J^cius = ^ ~ 8a/3 X 10-^. (172) 

Pc 

Assuming that af3 varies between 0.1 to 3, say, we get 

J^cius «i (s X 10-^ - 2.4 X 10-^) . (173) 



F{M, z) = erfc 



^/2aL{M,. 



erfc 



5c{l + z) 
V2ao{M) 
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Figure 6. (a) The Q, contributed by collapsed objects with mass greater than M plotted 
against a(M) at different values of z. The curves are for z = 0,1,2,3,4 and 5, from top 
to bottom. The constraint arising from cluster abundance at 2; = 0, quasar abundance at 
z = 2 — 4 and the abundance of damped Lyman-a systems at z = 2 — 3 are marked, (b) 
The M — a relation in a class of CDM-like models; 

[We shall concentrate on the top curve for z = 0, for the purpose of this 
example.] The fractional abundance given in (173), at z = 0, requires a 
a ~ (0.5 — 0.78) at the cluster scales. All wc need to determine now is 
whether a particular model has this range of a for cluster scales. Since 
this mass corresponds to a scale of about Sh'^Mpc, we conclude that the 
linearly extrapolated density contrast must be in the range cjl = (0.5 — 0.8) 
at R = 8h~^Mpc. This can act as a strong constraint on structure formation 
models. [The lower panel of figure (6) translates the bounds to a specific 
CDM model, parametrised by a shape parameter. This illustrates how any 
specific model can be compared with the bound in (171); for more details, 
see Padmanabhan, 1996 and references cited therein.] 

10. Scaling laws 

Before describing more sophisticated analytic approximations to gravita- 
tional clustering we shall briefly addres some simple scaling laws which can 
be obtained from our knowledge of linear evolution. These scaling laws are 
sufficiently powerful to allow reasonable predictions regarding the growth 
of structures in the universe and hence are useful for quick diagnostics of 
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a given model. We shall confine our attention to the scaling relations for 
a power-law spectrum for which \Sk\^ oc A;", and cr|^(-R) oc cx: 
j^-{n+3)/3_ |3ggj]^ \yy asking what restrictions can be put on the 

index n. 

The integrand defining in (95) behaves as near k = 0. [Note 

that Wk — 1 for small k in any window function]. Hence the finiteness 
of cj^ will require the condition n > —3. The behaviour of the integrand 
for large values of k depends on the window function W^- If we take the 
window function to be a Gaussian, then the convergence is ensured for all 
n. This might suggest that n can be made as large as one wants; that is, 
we can keep the power at small k (i.e., large wavelengths) to be as small as 
we desire. This result, however, is not quite true for the following reason: 
As the system evolves, small scale nonlinearities will develop in the system 
which can actually affect the large scales. If the large scales have too little 
power intrinsically (i.e. if n is large), then the long wavelength power will 
soon be dominated by the "tail" of the short wavelength power arising 
from the nonlinear clustering. This occurs because, in equation (36), the 
nonlinear terms and -Bk can dominate over AiTGpbS]^ at long wavelengths 
(as k ^ 0). Thus there will be an effective upper bound on n. 

The actual value of this upper-bound depends, to some extent, on the 
details of the small scale physics. It is, however, possible to argue that the 
natural value for this bound is n = 4. The argument runs as follows: Let us 
suppose that a large number of particles, each of mass m, are distributed 
carefully in space in such a way that there is very little power at large 
wavelengths. [That is, \6k\'^ oc A;" with n ^ 4 for small k]. As time goes on, 
the particles influence each other gravitationally and will start clustering. 
The density p(x, t) due to the particles in some region will be 

p(x,i) = ^m5[x-xi(t)], (174) 

i 

where Xj (t) is the position of the i-th particle at time t and the summation 
is over all the particles in some specified region. The density contrast in the 
Fourier space will be 

'^k(0 = ^E(exp[^k.x,(t)]-l) (175) 

i 

where N is the total number of particles in the region. For small enough 
|k|, we can expand the right hand side in a Taylor series obtaining 

S^it) = ik. 1 1 ^ X, (t) I - ^ 1 1 ^ xf (t) I + • • • . (176) 
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If the motion of the particles is such that the centre-of-mass of each of the 
subregions under consideration do not change, then J2 will vanish; under 
this (reasonable) condition, (5k (t) oc k'^ for small k. Note that this result 
follows, essentially, from the three assumptions: small-scale graininess of the 
system, conservation of mass and conservation of momentum. This will lead 
to a long wavelength tail with |(5jfcp oc /c^ which corresponds to n = 4. The 
corresponding power spectrum for gravitational potential Pip{k) oc /c^^jj^p 
is a constant. Thus, for all practical purposes, —3 < n < 4. The value n = 4 
corresponds to a\j{R) oc oc M~^/^. For comparison, note that purely 
Poisson fluctuations will correspond to {SM/M)'^ oc (1/M); i.e. (Tj^{R) oc 
oc with an index of n = 0. 
A more formal way of obtaining the A;^ tail is to solve equation (113) 
for long wavelengths; i.e. near k = (Padmanabhan, 1998). Writing (f)]^^ = 

+ (j)^^ + .... where (f)^^ = cf)^^ is the time independent gravitational 

(2) 

potential in the linear theory and (f)\^ is the next order correction, we get 
from (113), the equation 

(2) 

Writing (f)\^ = aCk one can determine Ck from the above equation. Plug- 
ging it back, we find the lowest order correction to be. 

Near k 2± 0, we have 

(2) ^ _Ja_ f^,:L 

oo 


(179) 

which is independent of k to the lowest order. Correspondingly the power 
spectrum for density Ps{k) oc k^Pip{k) oc k^ in this order. 

The generation of long wavelength k'^ tail is easily seen in simulations 
if one starts with a power spectrum that is sharply peaked in |k|. Figure 7 
shows the results of such a simulation (see Bagla and Padmanabhan, 1997) 
in which the y-axis is [A(fc)/a(t)]. In linear theory A oc a and this quantity 
should not change. The curves labelled by a = 0.12 to a = 20.0 shows the 
effects of nonlinear evolution, especially the development of A;^ tail. 



7,2 3 2 

8^ + 2^ 



5(k ■ p) 
A;2 
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Figure 7. The transfer of power to long wavelengths forming a fc* tail is illustrated 
using simulation results. Power is injected in the form of a narrow peak at L = 8 and 
the growth of power over and above the linear growth is shown in the figure. Note that 
the y — axis is (A/a) so that there will be no change of shape under linear evolution 
with A oc a. As time goes on a fc* tail is generated which itself evolves according to a 
nonlinear scaling relation discussed later on. 



Some more properties of the power spectra with different values of n can 
be obtained if the nonhnear effects are taken into account. We know that, 
in the matter-dominated phase, hnear perturbations grow as 5k{t) oc a(t) oc 
t^/^. Hence may assume that the perturbations 

at some scale R becomes nonlinear when aM{R) — 1. It follows that the 
time tR at which a scale R becomes nonlinear, satisfies the relation 

tR OC ii3(«+3)/4 ^ ^/(«+3)/4_ (-LgQ^ 

For n > —3, the timescale Ir is an increasing function of M; small scales 
become nonlinear at earlier times. The proper size L of the region which 

becomes nonlinear is 

L oc Ra{tR) oc oc i?(5+")/2 oc m(5+")/6. (181) 

Further, the objects which are formed at t = tR will have density p of the 
order of the background density p of the universe at tR. Since p oc f'^, we 
get 

p oc ^ i?-3(3+n)/2 ^ jy-(3+n)/2_ (^^82) 



48 



T. PADMANABHAN 



Combining (181) and (182) we get p oc L'l^ with 

/3=&tZ^. (183) 

In the nonlinear case, one may interpret the correlation function ^ as 
^(L) oc p{L); this would imply ^(x) oc x"!^ . ( We shall see later that such 
a behaviour is to be expected on more general grounds.) The gravitational 
potential due to these bodies is 

4) ~ Gp{L)L^ oc oc m(^-'^)/6. (184) 

The same scaling, of course, can be obtained from (j) oc (M/L). This re- 
sult shows that the binding energy of the structures increases with M for 
n < 1. In that case, the substructures will be rapidly erased as larger and 
larger structures become nonlinear. For n = 1, the gravitational potential 
is independent of the scale, and p oc L~^. 



11. Nonlinear scaling relations 

Given an initial density contrast, one can trivially obtain the density con- 
trast at any later epoch in the linear theory. If there is a procedure for 
relating the nonlinear density contrast and linear density contrast (even 
approximately) then one can make considerable progress in understanding 
nonlinear clustering. It is actually possible to make one such ansatz along 
the following lines. 

Let Vj.ci{a,x) denote the relative pair velocities of particles separated by 
a distance x, at an epoch a, averaged over the entire universe. This relative 
velocity is a measure of gravitational clustering at the scale x at the epoch 
a. Let h{a,x) = — [vj.f.\{a,x)/ax\ denote the ratio between the relative pair 
velocity and the Hubble velocity at the same scale. In the extreme nonlin- 
ear limit (^ ^ 1) , bound structures do not expand with Hubble flow. To 
maintain a stable structure, the relative pair velocity Urei (a, x) of parti- 
cles separated by x should balance the Hubble velocity Hr = dx; hence, 
Vj-ei = —dx or h (a, x) = 1. 

The behaviour of h (a, x) for ^ ^ 1 is more complicated and can be 
derived as follows: Let the peculiar velocity field be v (x) [wc shall suppress 
the a dependence since we will be working at constant a]. The mean relative 
velocity at a separation r = (x — y) is given by 

vrei(r) = ([v(x)-v(y)][l + 5(x)][l+J(y)]) 

- ([v(x)-v(y)]5(x)) + ([v(x)-v(y)]5(y)) (185) 

to lowest order, since S'^ term is higher order and (v (x) — v (y)) = 0. 
Denoting (v (x) — v (y)) by Vxy and writing x = y+r, the radial component 
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of relative velocity is 

'Vxy-r = j v(k) .1 



,jk.(r+y) _ ik.y 



(186) 



where v (k) is the Fourier transform of v (x) . This quantity is related to 5k 

by 

v(k) = i//a(^^^k. (187) 

(This equation is same as u = —Vip, used in (105), expressed in Fourier 
space). Using this in (186) and writing 6 (x) , 6 (y) in Fourier space, we find 
that 



v^y.r [6 (x) + (5 (y 
(fik f (fp 



iHa f [ ^ (H\ 5^,5* e'^^-'^yy 



gik.r _ ^ 



(188) 

We average this expression using {5)^5^) = (27r)'^ 6d (k — p) -P (fc) , to obtain 
Vrei-r = {vr,y-r[5{x)+5{y)]) 



d^k P{k) 



gik.r _ g-ik.r 



-2Ha 



J 



(Pk P(k) 
(27r)^ A;2 



(k.r) sin (k.r) , 



(189) 



From the symmetries in the problem, it is clear that Vrei(r) is in the direc- 
tion of r. So Vrei • r = Vre\r. The angular integrations are straightforward 
and give 



oo 

Ha f dk 

^rei = (vxy-r [5 (x) + 6 (y)]) = — 5- / —P (k) [kr cos kr — sin kr] . (190) 

rTT^ J k 



Using the expression (98) for ^ (r) this can be written as 

rv,,i{r) = ~(Har^)C (191) 
Dividing by r and noting that Hrpj-op = Har, we get 

(192) 



^ ^ ^'rel jr) ^ ^rel (r) ^ 2 - 



prop 



aHr 
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We get the important result that h{a, x) depends on (a, x) only through 
^(a, x) in the linear limit, while h = —lis the nonlinear limit. This suggests 
the ansatz that h depends on a and x only through some measure of the 
density contrast at the epoch a at the scale x. As a measure of the density 
contrast we shall use f(a,x) itself since the result in (192) clearly singles 
it out. In other words, we assume that h{a,x) = h[(,{a,x)]. Wc shall now 
obtain an equation connecting h and ^. By solving this equation, one can 
relate ^ and ^l. (Nityananda and Padmanabhan, 1994). The mean number 
of neighbours within a distance x of any given particle is 



Jo 



(193) 



when n is the comoving number density. Hence the conservation law for 
pairs implies 



dt ax^ dx 



[x'^{l+Ov]=0 



(194) 



where v{t,x) denotes the mean relative velocity of pairs at separation x 
and epoch t (We have dropped the subscript 'rel' for simplicity) . Using 



in (194), we get 
1 d 



d 



3x2 L 



(1+6] 



]_d_ 

dx 



ax- 



V d 
3 dx 



[^'(1 + 0] 



(195) 



(196) 



Integrating, we find: 



(197) 



[The integration would allow the addition of an arbitrary function of t on 
the right hand side. We have set this function to zero so as to reproduce 
the correct limiting behaviour] . It is now convenient to change the variables 
from t to a, thereby getting an equation for ^: 



-^^[x3(l+^-(a,x))] (198) 



ax J x^ dx 



or, defining h{a,x) = —{v/ax) 

d , d 
h- 



9 In a (? In X 



(1+0 = 3/1(1 + ^"). 



(199) 
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This equation shows that the behaviour of £,(a,x) is essentially decided by 
h, the dimensionless ratio between the mean relative velocity v and the 
Hubble velocity ax = {a/a)xprop, both evaluated at scale x. We shall now 
assume that 

h{x,a) = h[^{x,a)]. (200) 

This assumption, of course, is consistent with the extreme linear limit h = 
(2/3)^ and the extreme nonlinear limit h = 1. When h{x,a) = h[^{x,a)], it 
is possible to find a solution to (200) which reduces to the form ^ oc for 
^< 1 as follows: Let A = Ina, X = Inx and D{X,A) = (1 + 0- We define 
curves ("characteristics") in the X, A, D space which satisfy 



dX 



-h[D[X,A]] (201) 



i.e., the tangent to the curve at any point {X,A,D) is constrained by the 
value of h at that point. Along this curve, the left hand side of (199) is a 
total derivative allowing us to write it as 

This determines the variation of D along the curve. Integrating 

Squaring and determining the constant from the initial conditions at oq, in 
the linear regime 

... 1 1 f''^^ \ - 4 = (204) 



We now need to relate the scales x and I. Equation (201) can be written, 
using equation (202) as 

dX ^ 1 dD , , 

-TT = -h= TTT^-r-r (205 
dA 3DdA ^ ^ 

giving 

SX + hiD = ln[x^{l + ^)] = constant. (206) 
Using the initial condition in the linear regime, 



x^{l + ^)=l\ 



(207) 
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This shows that should be evaluated at I = It can be checked 

directly that (207 and (204) satisfy (199). The final result can, therefore be 
summarized by the equation (equivalent to (204) and (207)) 

/ 2 r^i"-'^) dp, \ - 
a(<..i)=exp(-/ S(rt(fTrtj^ i = .(l + f(a.x))V3. (208) 

Given the function ^(^), this relates and ^ or — equivalently — gives the 
mapping ^(a, x) = [/[^^(a, /)] between the nonlinear and linear correlation 
functions evaluated at different scales x and I. The lower limit of the integral 
is chosen to give In^ for small values of ^ on the linear regime. It may be 
mentioned that the equation (205) and its integral (207) are independent 
of the ansatz h{a,x) = h[^{a,x)]. 

The following points need to be stressed regarding this result: (i) Among 
all statistical indicators, it is only ^ which obeys a nonlinear scaling relation 
(NSR) of the form fNL(a, x) = U [^l(«, 0] • Attempts to write similar rela- 
tions for ^ or P{k) have no fundamental justification, (ii) The nonlocality of 
the relation represents the transfer of power in gravitational clustering and 
cannot be ignored — or approximated by a local relation between ^NL{a, x) 
and ^l{(i,x). 

Given the form of h{^), equation (208) determines the relation ^ = 
It is, however, easier to determine the form of U, directly along the 
following lines (Padmanabhan, 1996a): In the linear regime (^ <C 1,^l <C 1)) 
we clearly have U {^l) — ^l- To determine its form in the quasilinear regime, 
consider a region surrounding a density peak in the linear stage, around 
which we expect the clustering to take place. From the definition of ^ it 
follows that the density profile around this peak can be described by 

p{x) ^ Pbg[l + C{x)] (209) 

Hence the initial mean density contrast scales with the initial shell ra- 
dius / as (5j(/) oc ^l(/) in the initial epoch, when linear theory is valid. 
This shell will expand to a maximum radius of x^iax 
scale-invariant, radial collapse, models each shell may be approximated as 
contributing with a effective scale which is propotional to Xmax- Taking the 
final effective radius x as proportional to Xmax, the final mean correlation 
function will be 

M - 
^ql{x) oc p oc ^ oc ^p/^^^i^^s 0, a(0' (210) 

That is, the final correlation function in the quasilinear regime, ^qli at x 
is the cube of initial correlation function at I where oc x^^ oc x^£,ql{x). 
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Note that we did not assume that the initial power spectrum is a power 
law to get this result. In case the initial power spectrum is a power law, 
with x~^"'~^^\ then we immediately find that 

Cql OC x-^in+mn+i) (211) 

[If the correlation function in linear theory has the power law form S^l oc 
then the process described above changes the index from a to 3a/(l + a). 
We shall comment more about this aspect later]. For the power law case, 
the same result can be obtained by more explicit means. For example, in 
power law models the energy of spherical shell with mean density d(xi) oc 
x"^ will scale with its radius as £^ oc G5M{xi)/xi oc G5xf oc x^~^. Since 
M (X xf, it follows that the maximum radius reached by the shell scales as 
(M/E) oc . Taking the effective radius as a; = Xgff oc Xj"*" , the 
final density scales as 

In this quasilinear regime, ^ will scale like the density and we get ^ql oc 
^liQ index b can be related to n by assuming the the evolution 
starts at a moment when linear theory is valid. Since the gravitational 
potential energy [or the kinetic energy] scales as £' oc Xj in the linear 
theory, it follows that b = n + 3. This leads to the correlation function in 
the quasilinear regime, given by (211) . 

1{ Q = 1 and the initial spectrum is a power law, then there is no 
intrinsic scale in the problem. It follows that the evolution has to be self 
similar and ^ can only depend on the combination q = xa"^/^""'"^' . This 
allows to determine the a dependence of ^ql by substituting q for x in 
(211). We find 

^QL{a, X) oc o6/(n+4)^-3(n+3)/(n+4) (213) 

We know that, in the linear regime, ^ = oc a^. Equation (213) shows 
that, in the quasilinear regime, ^ = ^ql oc a^/("+^). Spectra with n < — 1 
grow faster than a^, spectra with n > — 1 grow slower than and n = — 1 
spectrum grows as a^. Direct algebra shows that 

^Qi(a,x)oc[a(a,Of (214) 

reconfirming the local dependence in a and nonlocal dependence in spatial 
coordinate. This result has no trace of original assumptions [spherical evo- 
lution, scale-invariant spectrum ....] left in it and hence once would strongly 
suspect that it will have far general validity. 
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Let us now proceed to the fully nonlinear regime. If we ignore the ef- 
fect of mergers, then it seems reasonable that virialised systems should 
maintain their densities and sizes in proper coordinates, i.e. the clustering 
should be "stable" . This would require the correlation function to have the 
form ^nl{cl, x) = a^F{ax). [The factor arising from the decrease in back- 
ground density] . From our previous analysis we expect this to be a function 
of ^L{a,l) where ^ x^^NL{a,x). Let us write this relation as 

^NL{a,x) = a^F{ax) = U^aJ)] (215) 

where U[z] is an unknown function of its argument which needs to be deter- 
mined. Since linear correlation function evolves as we know that we can 
write ^1,(0,/) = a^(5[?^] where Q is some known function of its argument. 
[We are using rather than I in defining this function just for future con- 
venience of notation]. In our case = x'^^NL{a, x) = {ax)'^F{ax) = r'^Fir) 
where we have changed variables from (a, x) to (a, r) with r = ax. Equation 
(215) now reads 

a^'Fir) = U[h{a, I)] = Uia'Qf]] = U[a'Q[r''F{r)]] (216) 

Consider this relation as a function of a at constant r. Clearly we need to 
satisfy [/"[cia^] = C2a^ where c\ and C2 are constants. Hence we must have 

U[z] oc z^/^ (217) 
Thus in the extreme nonlinear end we should have 

e7VL(a,x)oc [a(a,0]'/' (218) 

[Another way deriving this result is to note that if ^ = a'^F(ax), then 
h = 1. Integrating (208) with appropriate boundary condition leads to 
(218). ] Once again we did not need to invoke the assumption that the 
spectrum is a power law. //it is a power law, then we get, 

eiVL(a, x) cx aMa;-T; ^ = ^-^^ (219) 

(n + 5) 

This result is based on the assumption of ''stable clustering" and was orig- 
inally derived by Peebles (Peebles, 1980). It can be directly verified that 
the right hand side of this equation can be expressed in terms of q alone, 
as we would have expected. 

Putting all our results together, we find that the nonlinear mean cor- 
relation function can be expressed in terms of the linear mean correlation 
function by the relation: 



^(a,x) 



a(a,0_ (fora<l, f<l 

[ 14.14^^^(0, if'^ (for 5.85 < a, 200 < C) 



a(a, 0'^ (for 1< a < 5.85, 1< ^ < 200) (220) 
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The numerical coefficients have been determined by continuity arguments. 
We have assumed the hnear result to be valid upto ^ = 1 and the virialisa- 
tion to occur at ^ ~ 200 which is result arising from the spherical model. 
The exact values of the numerical coefficients can be obtained only from 
simulations. 

The true test of such a model, of course, is N-body simulations and 

remarkably enough, simulations are very well represented by relations of 
the above form. The simulation data for CDM, for example, is well fitted 
by (Padmanabhan etal., 1996): 



r^KO (for a < 1-2, e< 1-2) 



a (a, if (for 1< a < 5, 1< C < 125) (221) 

t 11.7a (a, i f^ (for 5 < a, 125 < 



which is fairly close to the theoretical prediction. [The fact that numerical 
simulations show a correlation between ^{a,x) and ^L{a,l) was originally 
pointed out by Hamilton et al., (1991) who, however, tried to give a mul- 
tiparameter fit to the data. This fit has somewhat obscured the simple 
physical interpretation of the result though has the virtue of being accu- 
rate for numerical work.] 

A comparison of (220) and (221) shows that the physical processes which 
operate at different scales are well represented by our model. In other words, 
the processes descibed in the quasilinear and nonlinear regimes for an indi- 
vidual lump still models the average behaviour of the universe in a statistical 
sense. It must be emphasized that the key point is the "flow of informa- 
tion" from / to X which is an exact result. Only when the results of the 
specific model are recast in terms of suitably chosen variables, we get a 
relation which is of general validity. It would have been, for example, incor- 
rect to use spherical model to obtain relation between linear and nonlinear 
densities at the same location or to model the function h. 

It may be noted that to obtain the result in the nonlinear regime, we 
needed to invoke the assumption of stable clustering which has not been de- 
duced from any fundamental considerations. In case mergers of structures 
are important, one would consider this assumption to be suspect (see Pad- 
manabhan et al., 1996). We can, however, generalise the above argument 
in the following manner: If the virialised systems have reached stationarity 
in the statistical sense, the function h — which is the ratio between two 
velocities — should reach some constant value. In that case, one can inte- 
grate (208) and obatin the result ^nl = a^^F{a^x) where h now denotes 
the asymptotic value. A similar argument will now show that 



eiVL(a,x)oc[a(a,Op/' 



(222) 
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in the general case. For the power law spectra, one would get 

fKx).«-; . = 5^|±|^ (223) 

Simulations are not accurate enough to fix the value of h; in particular, 
the asymptotic value of h could depend on n within the accuracy of the 
simulations. It may be possible to determine this dependence by modelling 
mergers in some simplified form. 

If /i = 1 asymptotically, the correlation function in the extreme nonlin- 
ear end depends on the linear index n. One may feel that physics at highly 
nonlinear end should be independent of the linear spectral index n. This 
will be the case if the asymptotic value of h satisfies the scaling 

h = (224) 
n + 3 

in the nonlinear end with some constant c. Only high resolution numerical 
simulations can test this conjecture that h{n + 3) = constant. 

It is possible to obtain similar relations between $,{a,x) and ^l(o, in 
two dimensions as well by repeating the above analysis (see Padmanabhan, 
1997). In 2-D the scaling relations turn out to be 

( i_L{a,l) (Linear) 
^(a,x) oc < ,f/^(a, /)^ (Quasi-linear) (225) 
ihlaJf'^ (Nonlinear) 

where h again denotes the asymptotic value. For power law spectrum the 
nonlinear correction function will iNL{a,x) = o?~^x~^ with 7 = 2(n -|- 
2)/(n + 4). 

If we generalize the concept of stable clustering to mean constancy 
of h in the nonlinear epoch, then the correlation function will behave as 
^Nhip-iX) = a^^Fia^x). In this case, if the spectrum is a power law then 
the nonlinear and linear indices are related to 

jM!i±2L (226) 
' 2 + /i(n + 2) ^ ' 

All the features discussed in the case of 3 dimensions arc present here as 
well. For example, if the asymptotic value of h scales with n such that h{n + 
2) = constant then the nonlinear index will be independent of the linear 
index. Figure (8) shows the results of numerical simulation in 2D, which 
suggests that h = 3/2 asymptotically (Bagla etal., 1998) We shall now 
consider some applications and further generalisations of these nonlinear 
scaling relations. 
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Figure 8. The comparison between theory and simulations in 2D. 



The ideas presented here can be generahsed in two obvious directions 
(see Munshi and Padmanabhan, 1997): (i) By considering peaks of different 
heights, drawn from an initial gaussian random field, and averaging over the 
probability distribution one can obtain a more precise NSR. (ii) By using a 
generalised ansatz for higher order correlation functions, one can attempt 
to compute the Sn parameters in the quasilinear and nonlinear regimes. 
We shall briefly comment on the results of these two generalisations. 

(i) The basic idea behind the model used to obtain the NSR can be de- 
scribed as follows: Consider the evolution of density perturbations starting 
from an initial configuration, which is taken to be a realisation of a Gaus- 
sian random field with variance a. A region with initial density contrast 5i 
will expand to a maximum radius Xf = Xi/5i and will contribute to the two- 
point correlation function an amount proportional to (xi/xj)^ = 6f. The 
initial density contrast within a randomly placed sphere of radius Xi will be 
ua{xi) with a probability proportional to exp(— 1/^/2). On the other hand, 
the initial density contrast within a sphere of radius Xi, centered around a 
peak in the density field will be proportional to the two-point correlation 
function and will be v^^(xi) with a probability proportional to exp(— 
It follows that the contribution from a typical region will scale as ^ oc ^^^^ 
while that from higher peaks will scale as ^ oc ^f. In the quasilinear phase, 
most dominant contribution arises from high peaks and we find the scaling 
to be ^QL oc ^f. The non-linear, virialized, regime is dominated by contri- 
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bution from several typical initial regions and has the scaling ^atl oc . 
This was essentially the result obtained above, except that we took v = \. 
To take into account the statistical fluctuations of the initial Gaussian field 
we can average over different v with a Gaussian probability distribution. 

Such an analysis leads to the following result. The relationship between 
^(a,x) and ^^(a, /) becomes 



C(a,x) = A[a(a,0]"'/';^ 



?,h 
2 



20F 



Zhja 



where 



6/i 



a = 



(227) 



(228) 



2 + h{n + 3) 

and A ^ 0.5 is the ratio between the final virialized radius and the radius at 
turn-around. In our model, h = 2 in the quasi-linear regime, and /i = 1 in 
the non-linear regime. However, the above result holds for any other value 
of h. Equation (227) shows that the scaling relations (220) acquire coeffi- 
cients which depend on the spectral index n when we average over peaks 
of different heights. (Mo etal., 1995; Munshi and Padmanabhan, 1997). 

(ii) In attempting to generalize our results to higher order correlation 
functions, it is important to keep the following aspect in mind. The Nth. or- 
der correlation function will involve TV — 1 different length scales. To make 
progress, one needs to assume that, although there are different length 
scales present in reduced n-point correlation function, all of them have to 
be roughly of the same order to give a significant contribution. If the corre- 
lation functions are described by a single scale, then a natural generalisation 
will be 

e^^(xf^-^V^'^^-'^ (229) 
Given such an ansatz for the N point correlation function, one can compute 
the Sn coefficients defined by the relation Sn = ^ straightfor- 

ward manner. We find that 



Sn = (47r 



,(JV-2)/2 



r(^) 



iV-l 



(230) 



where a is defined in equation (228). Given the function h{^), this equation 
allows one to compute (approximately) the value of Sn parameters in the 
quasi-linear and non-linear regimes. In our model h = 2 in the quasi-linear 
regime and ^ = 1 in the non-linear regime. The numerical values of Sjsf 
computed for different power spectra agrees reasonably well with simulation 
results. (For more details, see Munshi and Padmanabhan, 1997.) 
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12. NSR and halo profiles 

Now that we have a NSR giving ^(a, x) in terms of ^L(a, I) we can ask the 
question: How does the gravitational clustering proceed at highly nonlinear 
scales or, cquivalently, at any given scale at large a ? 

To begin with, it is easy to see that we must have v = —ax or h = 1 
for sufficiently large ^(a, x) if we assume that the evolution gets frozen in 
proper coordinates at highly nonlinear scales. Integrating equation (208) 
with h = 1, we get ^{a,x) = a^F{ax); this is the phenomenon we called 
"stable clustering". There are two points which need to be emphasised 
about stable clustering: 

(1) At present, there exists some evidence from simulations (see Pad- 
manabhan etal., 1996) that stable clustering does not occur in a = 1 
model. In a formal sense, numerical simulations cannot disprove [or even 
prove, strictly speaking] the occurrence of stable clustering, because of the 
finite dynamic range of any simulation. 

(2) . Theoretically speaking, the "naturalness"' of stable clustering is of- 
ten overstated. The usual argument is based on the assumption that at very 
small scales — corresponding to high nonlinearities — the structures are 
"expected to be" frozen at the proper coordinates. However, this argument 
does not take into account the fact that mergers are not negligible at any 
scale in an J7 = 1 universe. In fact, stable clustering is more likely to be 
valid in models with < 1 — a claim which seems to be again supported 
by simulations (see Padmanabhan etal., 1996). 

//stable clustering is valid, then the late time behaviour of ^(a, x) can- 
not be independent of initial conditions. In other words the two require- 
ments: (i) validity of stable clustering at highly nonlinear scales and (ii) 
the independence of late time behaviour from initial conditions, are mutu- 
ally exclusive. This is most easily seen for initial power spectra which are 
scale-free. If Pinik) oc /c" so that ^^(a, x) oc a^x~^'^^^\ then it is easy to 
show that ^(a, x) at nonlinear scales will vary as 



if stable clustering is true. Clearly, the power law index in the nonlinear 
regime "remembers" the initial index. The same result holds for more gen- 
eral initial conditions. 

What does this result imply for the profiles of individual halos? To an- 
swer this question, let us start with the simple assumption that the density 
field p{a, x) at late stages can be expressed as a superposition of several 
halos, each with some density profile; that is, we take 



_6_ 3(n+3) 

^a, x) OC a"+5a; "+5 ; 



(e » 200) 



(231) 




(232) 
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where the i-th halo is centered at x, and contributes an amount /(x — Xj, a) 
at the location [We can easily generalise this equation to the situation 
in which there arc halos with different properties, like core radius, mass etc 
by summing over the number density of objects with particular properties; 
we shall not bother to do this. At the other extreme, the exact description 
merely corresponds to taking the /'s to be Dirac delta functions. Hence 
there is no loss of generality in (232)]. The power spectrum for the density 
contrast, 5(a,x) = {p/ ph — 1), corresponding to the /9(a, x) in (232) can be 
expressed as 



P(k,a) oc (a=^|/(k,a)|)' 



exp — ik • Xj(a) 



(233) 



oc (a3|/(k,a)|)'Pcent(k,a) (234) 

where Pcent (k, a) denotes the power spectrum of the distribution of centers 
of the halos. 

If stable clustering is valid, then the density profiles of halos are frozen 
in proper coordinates and we will have /(x — Xj, a) = /(a (x — Xj)); hence 
the fouricr transform will have the form /(k, a) = /(k/a). On the 
other hand, the power spectrum at scales which participate in stable clus- 
tering must satisfy P(k, a) = P(k/a) [This is merely the requirement 
^(a, x) = a^F{ax) re-expressed in fourier space]. Prom equation (234) it 
follows that we must have Pcent(k, a) = Pccnt(k/a). We can however take 
Pcent = constant at sufficiently small scales. This is because we must neces- 
sarily have Pcent ~ constant, (by definition) for length scales smaller than 
typical halo size, when we are essentially probing the interior of a single 
halo at sufficiently small scales. We can relate the halo profile to the corre- 
lation function using (234). In particular, if the halo profile is a power law 
with / oc r~^, it follows that the ^(a,.x) scales as x~"' [see also McClelland 
and Silk, 1977] where 

7 = 2e - 3 (235) 

Now if the correlation function scales as ,x[^'^("+'^)/("^"^)l , then we see that 
the halo density profiles should be related to the initial power law index 
through the relation 

. = H!L±i) (236) 

So clearly, the halos of highly virialised systems still "remember" the initial 
power spectrum. 

Alternatively, without taking the help of the stable clustering hypoth- 
esis, one can try to "reason out"' the profiles of the individual halos and 
use it to obtain the scaling relation for correlation functions. One of the 
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favourite arguments used by cosmologists to obtain such a "reasonable" 
halo profile is based on spherical, scale invariant, collapse. It turns out that 
one can provide a scries of arguments, based on spherical collapse, to show 
that — under certain circumstances — the density profiles at the nonlinear 
end scale as The simplest variant of this argument runs as 

follows: If we start with an initial density profile which is r"", then scale 
invariant spherical collapse will lead to a profile which goes as r^'^ with 
P = 3a/(l + a) [see eg., Padmanabhan, 1996a]. Taking the intial slope as 
a = (n + 3)/2 will immediately give /? = 3(n + 3)/(n + 5). [Our definition of 
the stable clustering in the last section is based on the scaling of the corre- 
lation function and gave the slope of [— 3(ri + 3)/(n + 5)] for the correlation 
function. The spherical collapse gives the same slope for halo profiles.] In 
this case, when the halos have the slope of e = 3{n + 3)/(n + 5), then the 
correlation function should have slope 

Once again, the final state "remembers" the initial index n. 

Is this conclusion true ? Unfortunately, simulations do not have suffi- 
cient dynamic range to provide a clear answer but there are some claims 
that the halo profiles are "universal" and independent of initial conditions. 
The theoretical arguments given above are also far from rigourous (in spite 
of the popularity they seem to enjoy!). The argument for correlation func- 
tion to scale as [— 3(n + 3)/(n + 5)] is based on the assumption oi h = 1 
asymptotically, which may not be true. The argument, leading to density 
profiles scaling as , is based on scale invariant spherical col- 

lapse which does not do justice to nonradial motions. Just to illustrate the 
situations in which one may obtain final configurations which are indepen- 
dent of initial index n, we shall discuss two possibilities: 

(i) As a first example we will try to sec when the slope of the correlation 
function is universal and obtain the slope of halos in the nonlinear limit 
using our relation (235) . Such a situation can develop if we assume that h 
reaches a constant value asymptotically which is not necessarily unity. In 
that case, we get ^{a,x) = a^'^Fla^x] where h now denotes the constant 
asymptotic value of of the function. For an initial spectrum which is scale- 
free power law with index n, this result translates to 

^{a,x) (xa'^x"^ (238) 



where 7 is given by 



_ 3h{n + 3) 
^-2 + h{n + 3) ^^"^^^ 
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We now notice that one can obtain a 7 which is independent of initial power 
law index provided h satisfies the condition h{n + 3) = c, a constant. In 
this case, the nonlinear correlation function will be given by 



Note that we are now demanding the asymptotic value of h to explicitly 
depend on the initial conditions though the spatial dependence of ^(a, x) 
does not. In other words, the velocity distribution — which is related to 
h — still "remembers" the initial conditions. This is indirectly reflected in 
the fact that the growth of ^(a, x) — represented by a^'^/((^+'^)("'+^)) — does 
depend on the index n. 

We emphasize the fact that the velocity distribution remembers the 
initial condition because it is usual (in published literature) to ignore the 
memory in velocity and concentrate entirely on the correlation function. It 
is not clear to us [or we suppose to anyone else] whether it is possible to 
come up with a clustering scenario in which no physical feature remembers 
the initial conditions. This could probably occur when virialisation has run 
its full course but even then it is not clear whether the particles which 
evaporate from a given potential well (and form a uniform hot component) 
will forget all the initial conditions. 

As an example of the power of such a — seemingly simple — analysis 
note the following: Since c > 0, it follows that e > (3/2); invariant profiles 
with shallower indices (for e.g with e = 1) are not consistent with the 
evolution described above. 

(ii) For our second example, we shall make an ansatz for the halo profile 
and use it to determine the correlation function. We assume, based on small 
scale dynamics, that the density profiles of individual halos should resemble 
that of isothermal spheres, with e = 2, irrespective of initial conditions. 
Converting this halo profile to correlation function in the nonlinear regime 
is straightforward and is based on equation (235): If e = 2, we must have 
7 = 2e — 3 = 1 at small scales; that is ^{a,x) oc at the nonlinear regime. 
Note that this corresponds to the index at the nonlinear end, for which the 
growth rate is — same as in linear theory. We shall say more about 
such 'critical' indices later. [This growth however, is possible for initial 
power law spectra, only if e = 2, i.e h{n + 3) = 1 at very nonlinear scales. 
Testing the conjecture that h{n + 3) is a constant is probably a little easier 
than looking for invariant profiles in the simulations but the results are still 
uncertain] . 

The corresponding analysis for the intermediate regime, with 1 ^ ^{ajx) ^ 
is more involved. This is clearly seen in equation (234) which shows that 
the power spectrum [and hence the correlation function] depends both on 




(240) 
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the fourier transform of the halo profiles as well as the power spectrum 
of the distribution of halo centres. In general, both quantities will evolve 
with time and we cannot ignore the effect of -Pcciit(^) c^) a-^d relate P{k,a) 
to f{k,a). The density profile around a local maxima will scale approxi- 
mately as p oc ^ while the density profile around a randomly chosen point 
will scale as p oc [The relation 7 = 2e - 3 expresses the latter scal- 

ing of ^ oc p^]. There is, however, reason to believe that the intermediate 
regime (with 1 ^ ^ ^ 200) is dominated by the collapse of high peaks 
(see Padmanabhan, 1996a) . In that case, we expect the correlation func- 
tion and the density profile to have the same slope in the intermediate 
regime with ^{a,x) oc (1/x^). Remarkably enough, this corresponds to the 
'critical' index nc = — 1 for the intermediate regime for which the growth 
is proportional to a^. 

We thus see that if: (i) the individual halos are isothermal spheres with 
(1/a;^) profile and (ii) if ^ oc p in the intermediate regime and ^ oc p^ in 
the nonlinear regime, we end up with a correlation function which grows as 

at all scales. Such an evolution, of course, preserves the shape and is a 
good candidate for the late stage evolution of the clustering. 

While the above arguments are suggestive, they are far from conclusive. 
It is, however, clear from the above analysis and it is not easy to provide 
unique theoretical reasoning regarding the shapes of the halos. The situation 
gets more complicated if we include the fact that all halos will not all 
have the same mass, core radius etc and we have to modify our equations 
by integrating over the abundance of halos with a given value of mass, 
core radius etc. This brings in more ambiguities and depending on the 
assumptions we make for each of these components [c.g, abundance for 
halos of a particular mass could be based on Press-Schecter formalism], 
and the final results have no real significance. It is, therefore, better [and 
probably easier] to attack the question based on the evolution equation for 
the correlation function rather than from "physical" arguments for density 
profiles. 

13. Power transfer and critical indices 

Given a model for the evolution of the power spectra in the quasilinear and 
nonlinear regimes, one could generalise the questions raised in the last sec- 
tion and explore whether evolution of gravitational clustering possesses any 
universal charecteristics. For example one could ask whether a complicated 
initial power spectrum will be driven to any particular form of power spec- 
trum in the late stages of the evolution. This is a somewhat more general 
issue than, say, the invariance of halo profile. 

One suspects that such a possibility might arise because of the following 
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reason: We saw in the section 11 that [in the quasiUnear regime] spectra with 
n < — 1 grow faster than while spectra with n > — 1 grow slower than o? . 
This feature could drive the spectral index to n = rig ~ —1 in the quasilincar 
regime irrespective of the initial index. Similarly, the index in the nonlinear 
regime could be driven to n Ri — 2 during the late time evolution. So the 
spectral indices —1 and —2 are some kind of "fixed points" in the quasilinear 
and nonlinear regimes. Speculating along these lines, we would expect the 
gravitational clustering to lead to a "universal" profile which scales as 
at the nonlinear end changing over to x~'^ in the quasilinear regime. 

This effect can be understood better by studying the "effective" index 
for the power spectra at different stages of the evolution (see Bagla and 
Padmanabhan, 19977). To do this most effectively, let us define a local 
index for rate of clustering by 

na{a,x)= ' 241 

a mo 

which measures how fast ^(a, x) is growing. When ^(a, x) <C 1, then na = 2 
irrespective of the spatial variation of ^(a, a;) and the evolution preserves 
the shape of ^(a,x). However, as clustering develops, the growth rate will 
depend on the spatial variation of ^(a, x). Defining the effective spatial 
slope by 

- a, X + 3^ 242 
a m X 

one can rewrite the equation (199) as 

na = /i(77^ - (243) 

At any given scale of nonlinearity, decided by ^(a, x), there exists a critical 
spatial slope nc such that > 2 for < ric [implying rate of growth is 
faster than predicted by linear theory] and Ua < 2 for Ux > ric [with the 
rate of growth being slower than predicted by linear theory]. The critical 
index rif. is fixed by setting Ua = 2 in equation (243) at any instant. This 
requirement is established from the physically motivated desire to have a 
form of the two point correlation function that remains invariant under 
time evolution. Since the linear end of the two point correlation function 
scales as a^, the required invariance of form constrains the index Ua to be 2 
at all scales . The fact that ria > 2 for Hx < ric and Ha < 2 for Hx > ric will 
tend to "straighten out" correlation functions towards the critical slope. 
[We are assuming that ^(a,x) has a slope that is decreasing with scale, 
which is true for any physically interesting case]. Prom the NSR it is easy 
to see that in the range 1 ^ ^ ^ 200, the critical index is Wc ~ — 1 and 
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for 200 ^ ^, the critical index is nc ~ —2. This clearly suggests that the 
local effect of evolution is to drive the correlation function to have a shape 
with (l/x) behaviour at nonlinear regime and in the intermediate 

regime. Such a correlation function will have ~ 2 and hence will grow 
at a rate close to a^. 

The three panels of figure (9) illustrate features related to the existence 
of fixed points in a clearer manner. In the top panel we have plotted index 
of growth Ua = {dln^{a,x)/dliia)x as a function of ^ in the quasilinear 
regime obtained from the best fit for NSR based on simulations. Curves 
correspond to an input spectrum with index n = —2, —1, 1. The dashed 
horizontal line at na = 2 represents the linear growth rate. An index above 
this dashed horizontal line will represent a rate of growth faster than linear 
growth rate and the one below will represent a rate which is slower than 
the linear rate. It is clear that - in the quasilinear regime - the curve for 
n = — 1 closely follows the linear growth while n = —2 grows faster and 
n = 1 grows slower; so the critical index IS Tin -L . 

The second panel of figure 9 shows the effective index function of 

the index n of the original linear spectrum at different levels of nonlinearity 
labelled by ^ = 1, 5, 10, 50, 100. We see that in the quasilinear regime, Ua > 
2 for n < —1 and < 2 for n > —1. 

The lower panel of figure 9 shows the slope rix = —3 — {d\n^^/dlnx)a 
of ^ for different power law spectra. It is clear that crowds around 
nc ~ — 1 in the quasilinear regime. If perturbations grow by gravitaional 
instability, starting from an epoch at which ^initial <S 1 at all scales, then 
equation (243) with no > requires, that n^,, at any epoch, must satisfy 
the inequality 

Tlx < (3/0- (244) 

This bounding curve is shown by a dotted line in the figure. This powerful 
inequality shows that regions of strong nonlinearity [with ^ 1] should 
have effective index which is close to or less than zero. The index nc = — 1 
corresponds to the isothermal profile with ^{a, x) = a^x~'^ and has two 
interesting features to recommend it as a candidate for fixed point: 

(i) For n = — 1 spectra each logarithmic scale contributes the same 
amount of correlation potential energy. If the regime is modelled by scale 
invariant radial flows, then the kinetic energy will scale in the same way. It 
is conceivable that flow of power leads to such an equipartition state as a 
flxed point though it is difficult prove such a result in any generality. 

(ii) It was shown earlier that scale invariant spherical collapse will 
change the density profile x~^ with an index h to another profile with index 
36/(1 + 6). Such a mapping has a nontrivial fixed point for 6 = 2 corre- 
sponding to the isothermal profile and an index of power specrum n = — 1 
(see Padmanabhan, 1996a). 
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Figure 9. The top panel shows exponent of rate of growth of density fluctuations as a 
function of amplitude. We have plotted the rate of growth for three scale invariant spectra 
n = —2, —1, 1. The dashed horizontal line indicates the exponent for linear growth. For 
the range 1 < 5 < 100, the n = —1 spectrum grows as in linear theory; n < — 1 grows 
faster and n > — 1 grows slower. The second panel shows exponent of rate of growth as a 
function of linear index of the power spectrum for different values of ^ (1, 5, 10, 50, 100). 
Those arc represented by thick, dashed, dot-dashed, dotted and the dot-dot-dashed lines 
respectively. It is clear that spectra with riun < — 1 grow faster than the rate of growth 
in linear regime ajid nun > — 1 grow slower. The lower panel shows the evolution of index 
ria; = — 3 — {d\n(^/d\nx)a with ^. Indices vary from n = —2.5 to n = 4.0 in steps of 0.5. 
The tendency for nx to crowd around ric = — 1 is apparent_in the quasilinear regime. 
The dashed curve is a bounding curve for the index {ux < 3/^) if perturbations grow via 
gravitational instability. 
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These considerations also allow us to predict the nature of power trans- 
fer in gravitational clustering. Suppose that, initially, the power spectrum 
was sharply peaked at some scale /cq = 2tt/Lq and has a small width A/c. 
When the peak amplitude of the spectrum is far less than unity, the evolu- 
tion will be described by linear theory and there will be no flow of power 
to other scales. But once the peak approaches a value close to unity, power 
will be generated at other scales due to nonlinear couplings even though 
the amplitude of perturbations in these scales are less than unity. Mathe- 
matically, this can be understood from the evolution equation (36) for the 
density contrast — written in fourier space — as : 

4 + 2-4 = 47rGp(^k + Qk (245) 
a 

where 5]^{t) is the fourier transform of the density contrast, p is the back- 
ground density and Qk = — ^k is a nonlocal, nonlinear function which 
couples the mode k to all other modes k' . Coupling between different 
modes is significant in two cases. The obvious case is one with (5k > 1. A 
more interesting possibility arises for modes with no initial power [or ex- 
ponentially small power]. In this case nonlinear coupling provides the only 
driving terms, represented by Qk in equation (245). These generate power 
at the scale k through mode-coupling, provided power exists at some other 
scale. Note that the growth of power at the scale k will now he governed 
purely by nonlinear effects even though (5k <S 1. 

Physically, this arises along the following lines: If the initial spectrum is 
sharply peaked at some scale Lq, first structures to form are voids with a 
typical diameter Lq. Formation and fragmentation of sheets bounding the 
voids lead to generation of power at scales L < Lq. First bound structures 
will then form at the mass scale corresponding to Lq. In such a model, 
at L < Lq is nearly constant with an effective index of n ~ —3. Assuming 
we can use equation (220) with the local index in this case, we expect 
the power to grow very rapidly as compared to the linear rate of a^. [The 
rate of growth is for n = — 3 and for n = —2.5.] Different rate of 
growth for regions with different local index will lead to steepening of the 
power spectrum and an eventual slowing down of the rate of growth. In 
this process, which is the dominant one, the power transfer is mostly from 
large scales to small scales. [There is also a generation of the A;^ tail at large 
scales which wc have discussed earlier.] 

From our previous discussion, we would have expected such an evolution 
to lead to a "universal" power spectrum with some critical index nc ~ — 1 
for which the rate of growth is that of linear theory - viz., a^. In fact, the 
same results should hold even when there exists small scale power; recent 
numerical simulations dramatically confirm this prediction and show that 
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Figure 11. The growth of gravitational clustering towards a universal power spectrum 
P{k) oc A;-\ 
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- in the quasilinear regime, with 1 < 5 < 100 - power spectrum indeed 
has a universal slope (see figures 10, 11; for more details, see Bagla and 
Padmanabhan, 1997). 

The initial power spectrum for figure 10 was a Gaussian peaked at 
the scale ko = 2tt /Lq\Lq = 24 and having a spread Afc = 27r/128. The 
amplitude of the peak was chosen so that A;j„(fco = 27r/Lo,a = 0.25) = 1, 
where ^^{k) = k^P{k) /{2Tr'^) and P{k) is the power spectrum. Needless to 
say, the simulation starts while the peak of the Gaussian is in the linear 
regime (A(/co) ^ 1). The y-axis is A{k)/a, the power per logarithmic scale 
divided by the linear growth factor. This is plotted as a function of scale 
L = 2Tr/k for different values of scale factor a{t) and the curves are labeled 
by the value of a. As we have divided the power spectrum by its linear rate 
of growth, the change of shape of the spectrum occurs strictly because of 
non-linear mode coupling. It is clear from this figure that power at small 
scales grows rapidly and saturates to growth rate close to the linear rate 
[shown by crowding of curves] at later epochs. The effective index for the 
power spectrum approaches n = —1 within the accuracy of the simulations. 
Thus this figure clearly demonstrates the general features we expected from 
our understanding of scaling relations. 

Figure 11 compares power spectra of three different models at a late 
epoch. Model I was described in the last para; Model II had initial power 
concentrated in two narrow windows in /c-space. In addition to power 
around Lq = 24 as in model I, we added power at fei = 27r/Li; Li = 8 using 
a Gaussian with same width as that used in model I. Amplitude at Li was 
chosen five times higher than that at Lq = 24, thus A;j„(/ci, a = 0.05) = 1. 
Model III was similar to model II, with the small scale peak shifted to 
ki = 27r/Li; Li = 12. The amplitude of the small scale peak was the same 
as in Model II. At this epoch Auniko) = 4.5 and it is clear from this figure 
that the power spectra of these models are very similar to one another. 

There is another way of looking at this feature which is probably more 
useful. We recall that, in the study of finite gravitating systems made of 
point particles and interacting via newtonian gravity, isothermal spheres 
play an important role. They can be shown to be the local maxima of en- 
tropy [see Padmanabhan, 1990] and hence dynamical evolution drives the 
system towards an (1/x^) profile. Since one expects similar considerations 
to hold at small scales, during the late stages of evolution of the universe, 
we may hope that isothermal spheres with (1/a;^) profile may still play a 
role in the late stages of evolution of clustering in an expanding background. 
However, while converting the profile to correlation, we have to take note of 
the issues discussed earlier. In the intermediate regime, dominated by scale 
invariant radial collapse, the density will scale as the correlation function 
and we will have ^ cx (1/x^). On the other hand, in the nonlinear end, we 
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have the relation 7 = 2e — 3 which gives ^ oc for e = 2. Thus, if 

isothermal spheres are the generic contributors, then we expect the corre- 
lation function to vary as {l/x) and nonlinear scales, steepening to 
at intermediate scales. Further, since isothermal spheres are local maxima 
of entropy, a configuration like this should remain undistorted for a long 
duration. This argument suggests that a ^ which goes as (1/x) at small 
scales and at intermediate scales is likely to be a candidate for a 

pseudo-linear profile — that is configuration which grows approximately as 
at all scales. 

To go from the scalings in two limits to an actual profile, we can use some 
fitting function. By making the fitting function sufficiently complicated, we 
can make the pseudo-linear profile more exact. The simplest interpolation 
between the two limits is given by (Padmanabhan and Engineer, 1998) 



with L, B being constants. This approximate profile works reasonably well 
for the optimum value is B = 38.6. If we evolve this pseudo linear profile 
form = 1 to ~ 1000 using the NSR. and plot [^(a, a;)/a^] against x 
then the curves virtually fall on top of each other within about 10 per cent 
(see Padmanabhan and Engineer, 1998) This overlap of the curves show 
that the profile does grow approximately as a^. 

Finally, we will discuss a diff'erent way of thinking about pseudolinear 
profiles which may be useful. In studying the evolution of the density con- 
trast 5(a, x), it is conventional to expand in in term of the plane wave 
modes as 



In that case, the exact equation governing the evolution of (5(a,k) is given 



by the fact that in the linear regime we can ignore A and each of the modes 
evolve independently. For the same reason, this expansion is not of much 
value in the highly nonlinear regime. 

This prompts one to ask the question: Is it possible to choose some other 
set of basis functions (5(a,x), instead of exp ik • x, and expand (5(a,x) in 
the form 




(246) 




(247) 



k 



by 





(249) 



a 
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so that the nonUnear effects are minimised ? Here a stands for a set of pa- 
rameters describing the basis functions. This question is extremely difficult 
to answer, partly because it is ill-poscd. To make any progress, we have to 
first give meaning to the concept of "minimising the effects of nonlinear- 
ity". One possible approach we would like to suggest is the following: We 
know that when (5(a,x) ^ l,then ^(a,x) oc a F(x) for any arbitrary -F(x); 
that is all power spectra grow as in the linear regime. In the interme- 
diate and nonlinear regimes, no such general statement can be made. But 
it is conceivable that there exists certain special power spectra for which 
P(k, a) grows (at least approximately) as even in the nonlinear regime. 
For such a spectrum, the left hand side of (248) vanishes (approximately); 
hence the right hand side should also vanish. Clearly, such power spectra 
are affected least by nonlinear effects. Instead of looking for such a spe- 
cial P{k,a) wc can, cquivalcntly look for a particular form of S^{a,x) which 
evolves as closely to the linear theory as possible. Such correlation func- 
tions and corresponding power spectra [which are the pseudo-linear profiles] 
must be capable of capturing most of the essence of nonlinear dynamics. In 
this sense, wc can think of our pseudo-linear profiles as the basic building 
blocks of the nonlinear universe. The fact that the correlation function is 
closely related to isothermal spheres, indicates a connection between local 
gravitational dynamics and large scale gravitational clustering. 

14. Conclusion 

I tried to highlight in these lectures several aspects of gravitational cluster- 
ing which — I think — are important for understanding the basic physics. 
Some of the discussion points to obvious interrelationships with other branches 
of theoretical physics. For example, we saw that the power injected at any 
given scale cascades to other scales leading to a (nearly) universal power 
spectrum. This is remniscent of the fluid turbulence in which Kolmogorov 
spectrum arises as a (nearly) universal choice. Similarly, the existence of cer- 
tain configurations, which are least disturbed by the evolution [the "psuedo- 
linear profiles" , discussed in section 13] suggests similarities with the study 
of eddies in fluid mechanics, which possess a life of their own. Finally, the 
integral equation coupling the modes (113) promises to be an effective tool 
for analysing this problem. We are still far from having understood the dy- 
namics of this system from first principles and I hope these lectures serve 
the purpose of stimulating interest in this fascinating problem. 
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